Saturday, January 4, 2020

What Year Is It? VB6 Payload Crypter ! The year is 2020 !

What Year Is It? VB6 Payload Crypter
 JUNE 07, 2018
Last year, researchers identified new crimeware, Loki-Bot, which steals data and login credentials. Loki-Bot is generally distributed through malicious spam, and is difficult to identify without getting into the malware. Loki-Bot crimeware targets Windows, in contrast to the recent Android banking ransomware, Loki Bot. While the crimeware is not as prevalent in the wild, it has some unique and differentiating characteristics.


Loki-Bot’s crypter is especially interesting and unique because it utilizes Visual Basic 6.0 to load multiple stages of shellcode to deliver the Loki-Bot payload. We saw an interesting sample highlighted by @DissectMalware on Twitter and decided to take a closer look for ourselves. We’ll walk through Loki-Bot’s crypter functionality, the first and second stage shellcodes, the payload, and then provide some thoughts on stopping these kinds of attacks and what we can expect to see next.


Loki-Bot Crypter Functionality
The main Visual Basic compiled binary uses a raw pointer access technique to jump into the first stage shellcode which then calls the second-stage shellcode that is executed off of the stack. The first stage shellcode disables data execution prevention (DEP) to make the stack executable, and uses jmp esp to jump into the stage 2 shellcode to allow it to begin execution. Stage 2 then sets up persistence, decrypts the payload, and executes the payload contents by using process hollowing. In this case the payload is Loki-Bot, crimeware designed to steal private information from a system.  Loki-Bot’s functionality has already been covered in detail elsewhere, so we will instead focus on the mode of compromise and its anti-reversing engineering tactics.

DELIVERY CONTEXT
This particular malware sample began its life as an executable with an RTF exploit. Generally an RTF exploit is a specially crafted file that exploits vulnerabilities in the Rich Text Format parser of an application like Microsoft Word or Adobe Acrobat in order to gain code execution on the victim. The malware can then use this initial code execution to begin its exploit chain. Usually, crafted files are spread as attachments in phishing emails, which was the case in this sample. Once the malware gets code execution on the host it then downloads jazz.exe from the link below.
Source: hxxp://tpreiastephenville.com/jazz.exe
Sha256: a66f989e58ada2eff729ac2032ff71a159c521e7372373f4a1c1cf13f8ae2f0c

PE DESCRIPTION
The binary was compiled with VB6.0 professional/enterprise and so contains normal x86asm with a dependency on msvbvm60.dll. Stage 1 makes specific use of the visual basic runtime DLL to make dll calls to other libraries.
First Stage Shellcode
The first stage shellcode exists within the VB6 portion of the malware, which we’ll refer to as the crypter.  The first stage shellcode exists in the “rentegninger” sub module. The original sub module is then partially overwritten with obfuscated shellcode. The “Remanipulation8” public function is called from Load_Form(). This function manipulates the list of values of the Virtual Function Table returned by the “Me” reference to the form.
Set var_2 = Me

HOW TO GET THE SHELLCODE ENTRYPOINT
Stage 1 overwrites one of the variables in the compiled Visual Basic 6 code to point to an offset in the middle of the “rentegninger” submodule. The pointer is a hard coded integer that is calculated with division and a square root  as shown below.
var_17 = 0x1526C77
new_value = var_17 \ CLng(Sqr((25)))
new_value = 0x43AF4B
 
UTILIZING STRPTR TO ACCESS ADDRESSES
In the code example below the pointer to the Me value points to the beginning of the Virtual Function Table of the class that Me points to. In this case the Class is a Form. The offset 0x2B0 is actually the function “Show”. The pointer to the show function is overwritten by the entry point of the shellcode which is 0x43AF4B. Then you can easily call “Form.Show” and call into your shellcode. An example of this is located in the Appendix.
var_num1 = StrPtr(var_2) + 2B0h
ReplacePtr(var_num1, new_value)


This value is later called in the “Remanipulation8” function as call dword ptr [eax+2B0h]. It treats this call as a method of the Me object. DispCallFunc .
 

STRING-TO-STACK METHOD
The crypter uses a common trick to get the strings that are inline with the assembly onto the stack. When a call is made, the next address gets pushed onto the stack as the return destination address. A disassembler will try to disassemble the string even though it never gets executed.

USING SHELL_NOTIFYICON
The first stage shellcode uses the Shell_NotifyIcon in a non-standard way. It passes an address off the stack that does not resemble a proper PNOTIFYICONDATA struct. The Windows API still processes the events as if they were normal. As you can see below, the Shell_TrayWnd icon is junk data for this process. It is called twice where NIM_ADD and NIM_DELETE are used.



UTILIZING THE PEB LOADER AND DLLFUNCTIONCALL
The first stage shellcode uses a common technique among other Windows shellcode to get a reference to DllFunctionCall by utilizing the Process Environment Block (PEB). The PEB is a data structure provided to every running process, and can be used to gain information about that process such as environment variables, image base addresses and DLL imports. This shellcode contains a PEB loader routine that gets a reference to msvbvm60.dll  and then finds the offset of DllFunctionCall at 0x8D560CEC. Once it has the correct offset to DllFunctionCall, it can then use it to load Windows APIs so that it can make calls to them. More information on the PEB can be found here.
In a nutshell, the PEB is a linked list of offset values and the string names of the desired functions. You can linearly traverse the linked list, check for the desired function and save its offset if found.

DISABLING DATA EXECUTION PREVENTION (DEP)
The function ZwSetInformationProcess can be called with parameters -1 and 0x22 to turn off DEP for the process. DEP was originally intended to prevent programs from executing code on the stack, however DEP can also cause normal programs to crash without any notifications. Giving the program the ability to turn off DEP for itself allows the malware to avoid unknown crashes while also allowing malware authors to execute shellcode explicitly on the Stack. Microsoft Support provides additional details about DEP.



DECODING THE SHELLCODE
The stage 1 shellcode then utilizes JMP ESP to jump into that stack at the offset where the buffer for stage 2 shellcode was allocated. The stage 2 shellcode that was initially loaded onto the stack undergoes an initial pass of XOR decoding with an immediate value of 0x510473D1.


Second Stage Shellcode

SANDBOX EVASION


The stage 1 shellcode executes CPUID to detect if it’s being run in a virtual environment. If it is running in a virtual environment, it exits. If the stage 1 shellcode is not running in a virtual environment, then it continues execution normally. As detailed elsewhere, the malware sets the EAX register to 1, calls CPUID and then checks the 31st bit of the ECX register by applying a bitmask. If the 31st bit is 0 it knows it’s being run in a virtual environment.
Subsequently, the stage 1 shellcode calls the sleep function. Sleeping for prolonged durations of time is one evasion technique used by malware to subvert detection in sandboxed environments which are usually constrained by resource allocation to not run any given sample for longer than some established period of time, such as 30 seconds. Thus for the first 30 seconds of the program’s lifetime, it is benign and might fool some sandboxing environments.

ANTI-DEBUGGING USING NTYIELDEXECUTION
The ntdll function NtYieldExecution or its kernel32 equivalent SwitchToThread function allows the current thread to allocate the rest of the execution time, and allows the next scheduled thread to execute. If no threads are scheduled to execute or when the system is busy (and will not allow such a switch to occur), the ntdll NtYieldExecution() function returns a STATUS_NO_YIELD_PERFORMED (0x40000024) status code, which causes the kernel32 SwitchToThread function to return zero. When an application is being debugged, the act of single-stepping through the code causes debug events to occur and often results in no yield being allowed to occur. However, this is a hopelessly unreliable method of detecting the presence of a debugger because this method will also detect the presence of a thread that is running with high priority. An example of this code can be found here.

for (int i = 0; i < 0x20; i++)
{
         Sleep(0xf);
         if (NtYieldExecution() != STATUS_NO_YIELD_PERFORMED)
         iDebugged++;
}


CHECK FOR ADAPTERS AND WINDOWS
The stage 2 shellcode calls VirtualAllocEX to populate a new memory region with GetAdaptersInfo. It checks the offset +10Ch  of the struct for the Description. It then calls EnumWindows to check if the window has an empty string, most likely an attempt to detect execution within some sandbox.
PERSISTENCE AND PROCESS HOLLOWING
The stage 2 shellcode takes two routes. During the first route, the shellcode sets up persistent mechanisms with schtasks.exe. It then decrypts the payload with Xor and RC4 during the second route, creating a suspended process of itself and then hollows it out with the payload’s contents. Each route is explained below.
Route 1 (Persistence)
The stage 2 shellcode’s first route will acquire the hardcoded strings APPDATA=, TEMP=, and copied.exe in order to place a copy of itself in the %APPDATA% and %TEMP% locations as %AppData%\\Roaming\\copied.exe. Once the path is acquired, it will create a scheduled task to copy and run copied.exe using ShellExecuteA.
schtasks.exe" /Create /SC MINUTE /TN "Startup Key" /TR "%AppData%\\Roaming\\copied.exe"
schtasks.exe "/run /tn \"Startup Key\""

If that fails it will try again by adding "\" /RU SYSTEM"
schtasks.exe" /Create /SC MINUTE /TN "Startup Key" /TR "%AppData%\\Roaming\\copied.exe" /RU SYSTEM

It will also set the registry startup run key with:
schtasks.exe /Create /SC HOURLY /MO 12 /TN \"Startup Key\" /TR \"reg add \"HKLM\\Software\\Microsoft\\Windows\\CurrentVersion\\Run\" /v \"\\\"\"Startup Key\"\\\"\" /f /t REG_SZ /d \"\\\"\""

Route 2 (Process Hollowing)
The stage 2 shellcode’s second route focuses on decrypting the executable payload into memory and then hollowing out a child process to execute the payload. “Wee2" is the marker on the stack in the shellcode and in the file that denotes the beginning of  the copy operation. The format is Wee2<length of payload><key>.

Decrypting the Payload with Key
The key is stored in the shellcode after the marker and size in the shellcode. Its length is 0x100h, and is shown below.
00000000: 9944 4203 c046 b4f2 38dd 33ed 0281 473a  .DB..F..8.3...G:
00000010: 1c76 67d8 43bc d9c6 000f 58c2 c9f7 280e  .vg.C.....X...(.
00000020: 9fec 49ac 0bef bb56 8386 7d96 4c2a 4de3  ..I....V..}.L*M.
00000030: 221f 6e80 8e65 e02b 06b8 5f6a cf5c 72b7  ".n..e.+.._j.\r.
00000040: ea51 9354 1197 05ff 892e 843e 53d2 548b  .Q.T.......>S.T.
00000050: 6dc7 b829 940d e6d3 0d60 a913 d604 795f  m..).....`....y_
00000060: f0f9 9afd 183f 0ca7 d4d7 cee7 597b 9e34  .....?......Y{.4
00000070: 7370 bfd1 dfb6 317c 5709 b0bb 20ad c308  sp....1|W... ...
00000080: f7a2 e461 62e8 1250 da3b d54b a423 a5dc  ...ab..P.;.K.#..
00000090: be18 c536 e51a 3724 5eb1 fa20 2755 cab0  ...6..7$^.. 'U..
000000a0: 414a eb0a 6990 5df8 e1e4 dbf4 aacc ef41  AJ..i.]........A
000000b0: c4c1 10de ecc3 3ecd 645a 01c8 2dfe d015  ......>.dZ..-...
000000c0: 48f3 f1b2 b339 63a1 2b8c 269c f530 f6e9  H....9c.+.&..0..
000000d0: cb25 1687 366b 8875 af02 0771 78a6 1bbd  .%..6k.u...qx...
000000e0: 4e9b 3c5b bae1 ae49 3235 2c45 fbd9 fc92  N.<[...I25,E....
000000f0: 15ce 1d2f 3d14 8f1e b567 5219 7e4f 2166  .../=....gR.~O!f

The stage 2 shellcode uses the 0x100 sized byte key to xor the Loki-Bot payload of 0x34801 size then uses the same key to RC4 decrypt the product of the xor operation.

The shellcode then calls CreateProcessW to create a process of itself in suspended mode. It uses NTCreateFile, NTWriteVirtualMemory, NTUnMapViewofSection for process hollowing on the newly created process which contains the Loki-Bot payload. NtGetContextThread and NtSetContextThread and NtSetContextThread resume the process. The parent process terminates allowing the child to run as an orphan process. Because the Loki-Bot payload has already been reviewed by many researchers, we did not post the details about this malware.

Conclusion
There are a few key aspects to this crypter and its behaviour that make it fishy, including its crafty implementation of the VB6 runtime in shellcode, and use of anti-reverse engineering techniques and process hollowing. First, VB6 and the VB6 run time are rather old. While there are numerous binary distributions of software in the wild that were built with VB6 enterprise, it is still suspicious. Other suspicious activities include disabling its own DEP and checking if it’s being virtualized. Lastly, the crypter makes calls to the Windows API with malformed structs (i.e. the lack of an image for Shell_NotifyIcon). The combination of all of these suspicious activities could signal to a sensor like Endgame MalwareScoreTM that the program is up to no good, allowing us to stop it before the final execution occurs.  
As for the future, we are likely to see more samples using legacy run times and features. Judging from this sample, a performant Visual Basic 6 crypter has recently been distributed in the wild. It seems natural that in the future its capabilities will improve and the volume of distribution will increase with continued black market adoption.

Appendix: Stage 1 Code Replication
After careful analysis, we were able to reproduce the method of Stage 1 shellcode execution in VB6. The code below will only display an empty message box.

VB6 CODE
Private Declare Sub CopyMemory Lib "kernel32" Alias "RtlMoveMemory" _
    (ByRef lpvDest As Any, _
     ByRef lpvSrc As Any, _
     ByVal cbLength As Long)

Private Sub Form_Load()
    Me.Hide
    Dim lngRc As Long
    CopyMemory lngRc, Me, 4
    CopyMemory lngRc, ByVal lngRc, 4
    CopyMemory ByVal (lngRc + 688), 4202169, 4
    Me.Show
End Sub

SHELLCODE
E9 88 00 00 00 59 E9 8E 00 00 00 5A E8 0C 00 00 00 50 6A 00 6A 00 6A 00 6A 00
FF D0 C3 64 A1 30 00 00 00 8B 40 0C 8B 40 14 8B 00 8B 58 28 BE 4C 00 53 00 46
39 33 75 F1 81 7B 04 56 00 42 00 8B 70 10 56 8B 5E 3C 36 8B 34 24 01 DE 8B 5E
78 36 8B 04 24 01 D8 89 C6 83 C6 28 AD 85 C0 74 FB 03 04 24 BB 55 8B EC 83 39
18 75 EF 81 78 04 EC 0C 56 8D 75 E6 5B 31 DB 53 53 53 54 6A 00 81 04 24 00 00
04 00 52 51 54 FF D0 83 C4 1C C3 E8 73 FF FF FF 75 73 65 72 33 32 00 E8 6D FF
FF FF 4D 65 73 73 61 67 65 42 6F 78 41 00

ASSEMBLY
[SECTION .text]
global Start
Start:
    jmp     getdll
nextstring:
    pop ecx
    jmp     getstring
getapi:
    pop     edx ;put string in ebx
    call    loadapi
    push    eax ;get MessageBoxA
    push    0
    push    0
    push    0
    push    0
    call    eax
    retn
loadapi:
    mov     eax, [fs:0x30] ;Get the address of PEB
    mov     eax, [eax+0x0C] ;Get the address of PEB_LDR_DATA
    mov     eax, [eax+0x14] ;Get InMemoryOrderModuleList
loop1:
    mov     eax, [eax]
    mov     ebx, [eax+0x28]
    mov     esi, 0x53004C ;L S
    inc     esi ;M S for MSVBVM60.DLL
    cmp     [ebx], esi
    jnz     loop1
    cmp     dword [ebx+0x4], 0x420056 ;V B for MSVBVM60.DLL
    mov     esi, [eax+0x10]
    push    esi
    mov     ebx, [esi+0x3C]
    mov     esi, [ss:esp] ;<msvbvm60.Ordinal958>
    add     esi, ebx
    mov     ebx, [esi+0x78]
    mov     eax, [ss:esp] ;<msvbvm60.Ordinal958>
    add     eax, ebx
    mov     esi, eax
    add     esi, 0x28
loop2:
         lodsd
         test    eax, eax
         jz      short loop2
         add     eax, [esp]
         mov     ebx, 0x83EC8B55
         cmp     dword [eax], ebx
         jnz     short loop2
         cmp     dword [eax+0x4], 0x8D560CEC ;DllFunctionCall
         jnz     short loop2
         pop     ebx
         xor     ebx, ebx
         push    ebx
         push    ebx
         push    ebx
         push    esp
         push    0
         add     dword [esp], 0x40000
         push    edx ;MessageBoxA
         push    ecx ;user32
         push    esp
         call    eax ;DllFunctionCall
         add     esp, 0x1C
         retn
getdll:
    call nextstring
         db 0x75 ;u
         db 0x73 ;s
         db 0x65 ;e
         db 0x72 ;r
         db 0x33 ;3
         db 0x32 ;2
         db 0x00
getstring:
    call getapi
         db 0x4D ;M
         db 0x65 ;e
         db 0x73 ;s
         db 0x73 ;s
         db 0x61 ;a
         db 0x67 ;g
         db 0x65 ;e
         db 0x42 ;B
         db 0x6F ;o
         db 0x78 ;x
         db 0x41 ;A
         db 0x00


Source:
https://www.endgame.com/blog/technical-blog/what-year-it-vb6-payload-crypter






Thursday, January 2, 2020

Greatest common divisor (by Davide Chiappetta)

GCD (Greatest Common Divisor) of two numbers is the largest number that divides both of them. For example GCD of 20 and 28 is 4 and GCD of 98 and 56 is 14. In mathematics, the greatest common divisor (gcd) of two or more integers, which are not all zero, is the largest positive integer that divides each of the integers. For example, the gcd of 8 and 12 is 4.

Source:
https://davidechiappetta.wordpress.com/
http://davidechiappetta.blogspot.com/

Wednesday, December 25, 2019

Windows Information Provider (VB6 open source project)

Windows Information Provider Version 1.2.2 is an advanced app made by A_X_O (unfortunately I did not find his real name). It is a software that can teach you about everything inside the Windows 10 OS. This project is open source and can be modified, you can take parts of it as you need and you can ask A_X_O a little bit more about it (if he has the time to do so). You may wish to search for A_X_O on PSC/Internet, as he has published many advanced software as open source, made in the great VB6 programming language!

Download VB6 source code from PSC

Download VB6 from MediaFire

Download from ME

Final update: September 05th 2019
Now In Complete Applications Section

Drives: Drive, Label, File System, Serial...
Object Search: Object Name, DOS Path, Extension...
System: Description, Bootup State, Manufacturer...
Services: Service Name, ProcessID, State...
Accounts: Domain, Account, Password Required...
Software: Program Groups, Reliability Records, Installed Software...
Registry Search: Search Keys, Value, Data
Processes: Command Line, Handle, Parent Process ID, ProcessID...
TCP/UDP Tables: TCP/IP V4 & V6 Raw Data, UDP/IP V4 & V6 Raw Data
Peripherals: Keyboard, Pointing Device, Sound Card...
Recycle Bin Info: Item Count, Size
Processor: Current Clock Speed, Manufacturer...
Memory Info: BankLabel, Capacity, Manufacturer...
Security Settings: Anti-Malware, Anti-Spyware, Anti-Virus...
WMI Folder Listings: Details of Directories from WMI
WMI File Listings: Details of Files from WMI
View Hidden Object Status: View and Set Hidden Object Status
Common Folders: Quick Access to Common Directories
Google Chrome Extensions:Profile, Extension ID, Name, Version, Permissions, Installed, Key

Windows Information Provider Version 1.2.2













Sunday, December 8, 2019

Experimenting with Markov Chains: Prediction and Simulation

For those interested in experimenting with Markov Chains, we found a new type of app implemented in four or five different programming languages, including VB6. The applications below represent a prediction and a simulator for different experiments related to Markov Chains, both developed by Dr. Paul A. Gagniuc for a math book (Gagniuc, Paul A. (2017). Markov Chains: From Theory to Implementation and Experimentation. USA, NJ: John Wiley & Sons. pp. 1–235. ISBN 978-1-119-38755-8.). These apps are available as binary files and source code.


The prediction app:


Markov Chains - Prediction Framework

The simulation app:

Markov Chains - Simulation Framework

A weather prediction application with animation:


Download from me
Download from Wiley

Markov Chains - The Weather


Sources:

1. https://www.wiley.com/en-gb/Markov+Chains%3A+From+Theory+to+Implementation+and+Experimentation-p-9781119387558

2. https://transition-matrix.blogspot.com/




Discrete Probability Detector algorithm

'###################################################################################################
'# John Wiley & Sons, Inc.                                                                         #
'#                                                                                                 #
'# Book:   Markov Chains: From Theory To Implementation And Experimentation                        #
'# Author: Dr. Paul Gagniuc                                                                        #
'# Data:   01/05/2017                                                                              #
'#                                                                                                 #
'# Title:                                                                                          #
'# Discrete Probability Detector                                                                   #
'#                                                                                                 #
'# Short Description:                                                                              #
'# The purpose of this algorithm is to convert any string into a probability matrix.               #
'#_______________________                                                                          #
'# Detailed description: \                                                                         #
'# This algorithm is an advanced variation of the "ExtractProb" function from the book. The        #
'# main difference between "ExtractProb" function and the DPD algorithm is the automatic           #
'# identification of states.                                                                       #
'#_________________________________________________________________________________________________#
'# Initially, the states are identified in the first phase. Each new letter found in "S" is        #
'# appended to the string forming in variable "a". Thus, variable "a" gradually increases until    #
'# all types of letters from "S" are identified.                                                   #
'#_________________________________________________________________________________________________#
'# In the second phase the elements of matrix "m" are filled with zero values for later            #
'# purposes. Also, in the second phase the first column of matrix "e" is filled with letters       #
'# found in variable "a", and the second column of matrix "e" is filled with zero values for       #
'# later use.                                                                                      #
'#_________________________________________________________________________________________________#
'# In the third phase, the transitions between letters of "S" are counted and stored in            #
'# matrix "m".                                                                                     #
'# The strategy in this particular case is to fill matrix "m" with transition counts before the    #
'# last letter in "S" is reached. In this case, the first column of matrix "e" already contains    #
'# the letters from variable "a". The two components of vector "l" contain the "i" and "i+1"       #
'# letters from "S". The count of individual transitions between letters is made by a comparison   #
'# between vector "l" and the elements from the first column of matrix "e". The number of rows     #
'# in matrix "m" and matrix "e" is the same, namely "d". Therefore, an extra loop can be avoided   #
'# by mapping matrix "m" through a coordinate system. For instance, if the letter from position    #
'# "i" in "S" stored in "l(0)" and the letter from "j" row in matrix "e" (e(j,0)) are the same     #
'# then variable "r = j". Likewise, if the letter "i+1" stored in l(1) and the letter from e(j,0)  #
'# are the same then variable "c = j". Variable "r" represents the rows of matrix "m", whereas     #
'# variable "c" represents the columns of matrix "m" (m(r, c)).                                    #
'# Thus, at each step through "S", an element of matrix "m" is always incremented according to     #
'# the coordinates received from "r" and "c". This "coordinate" approach greatly increases the     #
'# processing speed of the algorithm. The number of loops = (k-1)*d, where "d" represents the      #
'# number of states (or letter types), and "k" is the number of letters in "S". When the letter    #
'# stored in "l(0)" and the letter from "j" row in matrix "e" are the same, the second column of   #
'# matrix "e" is also incremented. The second column of matrix "e" stores the number of            #
'# appearances for each type of letter in "S".                                                     #
'#_________________________________________________________________________________________________#
'# In the fourth phase, the counts from matrix "m" elements are divided by the counts from the     #
'# second column of matrix "e". The result of this division is stored in the same position in      #
'# matrix "m", and represents a transition probability.                                            #
'#_________________________________________________________________________________________________#
'#                                                                                                 #
'# Special considerations:                                                                         #
'#                                                                                                 #
'# If a state at the end of "S" (ie HAHAAAHQ) does not occur in the rest of "S" then matrix "m"    #
'# will contain a row with all elements on zero. Since it is at the end of "S", the letter does    #
'# not make a transition to anything. If a state from the beginning of "S" (ie. QHAHAAAH) does     #
'# not occur in the rest of "S" then matrix "m" will contain a column with all elements on zero.   #
'# Since the first letter it is only seen at the beginning of "S", no other letter makes a         #
'# transition to it.                                                                               #
'#                                                                                                 #
'# The meaning of variables:                                                                       #
'#      _____________________________________________________________________________________      #
'# S = |The string that is being analyzed.                                                  _|     #
'#      _____________________________________________________________________________________      #
'# q = |It is a flag variable with initial value of 1. The value of q becomes zero only if a |     #
'#     |letter x in the "S" string coresponds with a letter y in the "a" string.            _|     #
'#      _____________________________________________________________________________________      #
'# a = |The variable that holds the letters representing the states. The variable gradually  |     #
'#     |increases in length as the "S" string is analyzed. At each loop, a new letter is     |     #
'#     |added to variable "a" only if the value of q becomes zero. Thus, at the end of the   |     #
'#     |first loop the number of letters in the variable is equal to the total number        |     #
'#     |of states.                                                                          _|     #
'#      _____________________________________________________________________________________      #
'# d = |Represents the total number of states and is the length of "a" variable.            _|     #
'#      _____________________________________________________________________________________      #
'# m = |The main probability matrix which the function produces.                            _|     #
'#      _____________________________________________________________________________________      #
'# k = |Represents the length of the "S" string.                                            _|     #
'#      _____________________________________________________________________________________      #
'# e = |It is a matrix with two columns, namely column 0 and 1. Column 0 stores all the      |     #
'#     |letters found in "a". Column 1 stores the number of appearances for each type        |     #
'#     |of letter in "S".                                                                   _|     #
'#      _____________________________________________________________________________________      #
'# l = |Is a vector with two components. Vector "l" contains the "i" and "i+1" letters       |     #
'#     |from "S".                                                                           _|     #
'#                                                                                                 #
'###################################################################################################

Private Sub MakeMatrix_Click()
    Discrete_Probability_Detector (InP.Text)
End Sub

Function Discrete_Probability_Detector(ByVal S As String)

    Dim e() As String
    Dim m() As String
    Dim l(0 To 1) As String
    
    k = Len(S)
    w = 1
    
    For i = 1 To k
        q = 1
        For j = 0 To Len(a)
            x = Mid(S, i, 1)
            y = Mid(a, j + 1, 1)
            If x = y Then q = 0
        Next j
        If q = 1 Then a = a & x
    Next i

    d = Len(a)

    ReDim e(w To d, 0 To 1) As String
    ReDim m(w To d, w To d) As String
    
    For i = w To d
        For j = w To d
          m(i, j) = 0
          If j = w Then
            e(i, 0) = Mid(a, i, 1)
            e(i, 1) = 0
          End If
        Next j
    Next i

    For i = 1 To k - 1
        l(0) = Mid(S, i, 1)
        l(1) = Mid(S, i + 1, 1)
        For j = w To d
            If l(0) = e(j, 0) Then
               e(j, 1) = Val(e(j, 1)) + 1
               r = j
            End If
            If l(1) = e(j, 0) Then c = j
        Next j
        m(r, c) = Val(m(r, c)) + 1
    Next i
    
    tmp = "S=" & S & vbCrLf & vbCrLf
    tmp = tmp & "The algorithm detected a total of " & (d - w + 1) & " states." & vbCrLf & vbCrLf
    tmp = tmp & MatrixPaint(w, d, m, a, "(C)", "Count:") & vbCrLf
    
    For i = w To d
        For j = w To d
            If Val(e(i, 1)) > 0 Then
            m(i, j) = Val(m(i, j)) / Val(e(i, 1))
            End If
        Next j
    Next i

    tmp = tmp & MatrixPaint(w, d, m, a, "(P)", "Transition matrix M:")

    OutPut.Text = tmp

End Function


Function MatrixPaint(w, d, ByVal m As Variant, a, n, ByVal msg As String) As String
    
    Dim e() As String
    ReDim e(1 To d) As String
    
    d = Len(a)
    q = "|     "
    h = "|_____|"
    l = vbCrLf
    
    For i = (w - 1) To d
        If i = (w - 1) Then t = t & l & "."
        t = t & "_____."
        If i = d Then t = t & l & "| " & n & " |  "
    Next i

    For i = w To d
        e(i) = Mid(a, i, 1)
        t = t & e(i) & "  |  "
        h = h & "_____|"
    Next i
    
    t = t & l & h & l
    
    For i = w To d
        For j = w To d
            v = Round(m(i, j), 2)
            u = Mid(q, 1, Len(q) - Len(v))
            If j = d Then o = "|" Else o = ""
            For b = w To d
                If j = w And i = b Then
                    t = t & "|  " & e(i) & "  "
                End If
            Next b
            t = t & u & v & o
        Next j
    t = t & l & h & l
    Next i
    
    MatrixPaint = msg & " M[" & Val(d - w + 1) & "," & Val(d - w + 1) & "]" & l & t & l
    
End Function


Private Sub Form_Resize()
    If DPD.ScaleWidth > 0 Then
        OutPut.Width = DPD.ScaleWidth - OutPut.Left - 10
        OutPut.Height = DPD.ScaleHeight - OutPut.Top - 10
    End If
End Sub

Source: Gagniuc, Paul A. (2017). Markov Chains: From Theory to Implementation and Experimentation. USA, NJ: John Wiley & Sons. pp. 1–235. ISBN 978-1-119-38755-8.


Other sources:

1. https://sites.google.com/view/text-to-transition-matrix/

2. https://www.statisticsviews.com/details/feature/10822963/Markov-Chains-From-Theory-to-Implementation-and-Experimentation-An-interview-wit.html

3. https://www.wiley.com/en-gb/Markov+Chains%3A+From+Theory+to+Implementation+and+Experimentation-p-9781119387558

4. https://transition-matrix.blogspot.com/

Fast kernel

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Copyright (c) 2009, Sergey Bochkanov (ALGLIB project).
'
'>>> SOURCE LICENSE >>>
'This program is free software; you can redistribute it and/or modify
'it under the terms of the GNU General Public License as published by
'the Free Software Foundation (www.fsf.org); either version 2 of the
'License, or (at your option) any later version.
'
'This program is distributed in the hope that it will be useful,
'but WITHOUT ANY WARRANTY; without even the implied warranty of
'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
'GNU General Public License for more details.
'
'A copy of the GNU General Public License is available at
'http://www.fsf.org/licensing/licenses
'
'>>> END OF LICENSE >>>
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function CMatrixRank1F(ByVal M As Long, _
         ByVal N As Long, _
         ByRef A() As Complex, _
         ByVal IA As Long, _
         ByVal JA As Long, _
         ByRef U() As Complex, _
         ByVal IU As Long, _
         ByRef V() As Complex, _
         ByVal IV As Long) As Boolean
    Dim Result As Boolean

    Result = False

    CMatrixRank1F = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function RMatrixRank1F(ByVal M As Long, _
         ByVal N As Long, _
         ByRef A() As Double, _
         ByVal IA As Long, _
         ByVal JA As Long, _
         ByRef U() As Double, _
         ByVal IU As Long, _
         ByRef V() As Double, _
         ByVal IV As Long) As Boolean
    Dim Result As Boolean

    Result = False

    RMatrixRank1F = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function CMatrixMVF(ByVal M As Long, _
         ByVal N As Long, _
         ByRef A() As Complex, _
         ByVal IA As Long, _
         ByVal JA As Long, _
         ByVal OpA As Long, _
         ByRef X() As Complex, _
         ByVal IX As Long, _
         ByRef Y() As Complex, _
         ByVal IY As Long) As Boolean
    Dim Result As Boolean

    Result = False

    CMatrixMVF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function RMatrixMVF(ByVal M As Long, _
         ByVal N As Long, _
         ByRef A() As Double, _
         ByVal IA As Long, _
         ByVal JA As Long, _
         ByVal OpA As Long, _
         ByRef X() As Double, _
         ByVal IX As Long, _
         ByRef Y() As Double, _
         ByVal IY As Long) As Boolean
    Dim Result As Boolean

    Result = False

    RMatrixMVF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function CMatrixRightTRSMF(ByVal M As Long, _
         ByVal N As Long, _
         ByRef A() As Complex, _
         ByVal I1 As Long, _
         ByVal J1 As Long, _
         ByVal IsUpper As Boolean, _
         ByVal IsUnit As Boolean, _
         ByVal OpType As Long, _
         ByRef X() As Complex, _
         ByVal I2 As Long, _
         ByVal J2 As Long) As Boolean
    Dim Result As Boolean

    Result = False

    CMatrixRightTRSMF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function CMatrixLeftTRSMF(ByVal M As Long, _
         ByVal N As Long, _
         ByRef A() As Complex, _
         ByVal I1 As Long, _
         ByVal J1 As Long, _
         ByVal IsUpper As Boolean, _
         ByVal IsUnit As Boolean, _
         ByVal OpType As Long, _
         ByRef X() As Complex, _
         ByVal I2 As Long, _
         ByVal J2 As Long) As Boolean
    Dim Result As Boolean

    Result = False

    CMatrixLeftTRSMF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function RMatrixRightTRSMF(ByVal M As Long, _
         ByVal N As Long, _
         ByRef A() As Double, _
         ByVal I1 As Long, _
         ByVal J1 As Long, _
         ByVal IsUpper As Boolean, _
         ByVal IsUnit As Boolean, _
         ByVal OpType As Long, _
         ByRef X() As Double, _
         ByVal I2 As Long, _
         ByVal J2 As Long) As Boolean
    Dim Result As Boolean

    Result = False

    RMatrixRightTRSMF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function RMatrixLeftTRSMF(ByVal M As Long, _
         ByVal N As Long, _
         ByRef A() As Double, _
         ByVal I1 As Long, _
         ByVal J1 As Long, _
         ByVal IsUpper As Boolean, _
         ByVal IsUnit As Boolean, _
         ByVal OpType As Long, _
         ByRef X() As Double, _
         ByVal I2 As Long, _
         ByVal J2 As Long) As Boolean
    Dim Result As Boolean

    Result = False

    RMatrixLeftTRSMF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function CMatrixSYRKF(ByVal N As Long, _
         ByVal K As Long, _
         ByVal Alpha As Double, _
         ByRef A() As Complex, _
         ByVal IA As Long, _
         ByVal JA As Long, _
         ByVal OpTypeA As Long, _
         ByVal Beta As Double, _
         ByRef C() As Complex, _
         ByVal IC As Long, _
         ByVal JC As Long, _
         ByVal IsUpper As Boolean) As Boolean
    Dim Result As Boolean

    Result = False

    CMatrixSYRKF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function RMatrixSYRKF(ByVal N As Long, _
         ByVal K As Long, _
         ByVal Alpha As Double, _
         ByRef A() As Double, _
         ByVal IA As Long, _
         ByVal JA As Long, _
         ByVal OpTypeA As Long, _
         ByVal Beta As Double, _
         ByRef C() As Double, _
         ByVal IC As Long, _
         ByVal JC As Long, _
         ByVal IsUpper As Boolean) As Boolean
    Dim Result As Boolean

    Result = False

    RMatrixSYRKF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function RMatrixGEMMF(ByVal M As Long, _
         ByVal N As Long, _
         ByVal K As Long, _
         ByVal Alpha As Double, _
         ByRef A() As Double, _
         ByVal IA As Long, _
         ByVal JA As Long, _
         ByVal OpTypeA As Long, _
         ByRef B() As Double, _
         ByVal IB As Long, _
         ByVal JB As Long, _
         ByVal OpTypeB As Long, _
         ByVal Beta As Double, _
         ByRef C() As Double, _
         ByVal IC As Long, _
         ByVal JC As Long) As Boolean
    Dim Result As Boolean

    Result = False

    RMatrixGEMMF = Result
End Function


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Fast kernel
'
'  -- ALGLIB routine --
'     19.01.2010
'     Bochkanov Sergey
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function CMatrixGEMMF(ByVal M As Long, _
         ByVal N As Long, _
         ByVal K As Long, _
         ByRef Alpha_ As Complex, _
         ByRef A() As Complex, _
         ByVal IA As Long, _
         ByVal JA As Long, _
         ByVal OpTypeA As Long, _
         ByRef B() As Complex, _
         ByVal IB As Long, _
         ByVal JB As Long, _
         ByVal OpTypeB As Long, _
         ByRef Beta_ As Complex, _
         ByRef C() As Complex, _
         ByVal IC As Long, _
         ByVal JC As Long) As Boolean
    Dim Result As Boolean
    Dim Alpha As Complex
    Dim Beta As Complex
    Alpha = Alpha_
    Beta = Beta_

    Result = False

    CMatrixGEMMF = Result
End Function

Monday, July 23, 2018

Solution of the differential equation

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Cephes Math Library Release 2.8:  June, 2000
'Copyright by Stephen L. Moshier
'
'Contributors:
'    * Sergey Bochkanov (ALGLIB project). Translation from C to
'      pseudocode.
'
'See subroutines comments for additional copyrights.
'
'>>> SOURCE LICENSE >>>
'This program is free software; you can redistribute it and/or modify
'it under the terms of the GNU General Public License as published by
'the Free Software Foundation (www.fsf.org); either version 2 of the
'License, or (at your option) any later version.
'
'This program is distributed in the hope that it will be useful,
'but WITHOUT ANY WARRANTY; without even the implied warranty of
'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
'GNU General Public License for more details.
'
'A copy of the GNU General Public License is available at
'http://www.fsf.org/licensing/licenses
'
'>>> END OF LICENSE >>>
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Airy function
'
'Solution of the differential equation
'
'y"(x) = xy.
'
'The function returns the two independent solutions Ai, Bi
'and their first derivatives Ai'(x), Bi'(x).
'
'Evaluation is by power series summation for small x,
'by rational minimax approximations for large x.
'
'
'
'ACCURACY:
'Error criterion is absolute when function <= 1, relative
'when function > 1, except * denotes relative error criterion.
'For large negative x, the absolute error increases as x^1.5.
'For large positive x, the relative error increases as x^1.5.
'
'Arithmetic  domain   function  # trials      peak         rms
'IEEE        -10, 0     Ai        10000       1.6e-15     2.7e-16
'IEEE          0, 10    Ai        10000       2.3e-14*    1.8e-15*
'IEEE        -10, 0     Ai'       10000       4.6e-15     7.6e-16
'IEEE          0, 10    Ai'       10000       1.8e-14*    1.5e-15*
'IEEE        -10, 10    Bi        30000       4.2e-15     5.3e-16
'IEEE        -10, 10    Bi'       30000       4.9e-15     7.3e-16
'
'Cephes Math Library Release 2.8:  June, 2000
'Copyright 1984, 1987, 1989, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Sub Airy(ByVal x As Double, _
         ByRef Ai As Double, _
         ByRef Aip As Double, _
         ByRef Bi As Double, _
         ByRef Bip As Double)
    Dim z As Double
    Dim zz As Double
    Dim t As Double
    Dim f As Double
    Dim g As Double
    Dim uf As Double
    Dim ug As Double
    Dim k As Double
    Dim zeta As Double
    Dim theta As Double
    Dim domflg As Long
    Dim c1 As Double
    Dim c2 As Double
    Dim sqrt3 As Double
    Dim sqpii As Double
    Dim AFN As Double
    Dim AFD As Double
    Dim AGN As Double
    Dim AGD As Double
    Dim APFN As Double
    Dim APFD As Double
    Dim APGN As Double
    Dim APGD As Double
    Dim AN As Double
    Dim AD As Double
    Dim APN As Double
    Dim APD As Double
    Dim BN16 As Double
    Dim BD16 As Double
    Dim BPPN As Double
    Dim BPPD As Double

    sqpii = 0.564189583547756
    c1 = 0.355028053887817
    c2 = 0.258819403792807
    sqrt3 = 1.73205080756888
    domflg = 0#
    If x > 25.77 Then
        Ai = 0#
        Aip = 0#
        Bi = MaxRealNumber
        Bip = MaxRealNumber
        Exit Sub
    End If
    If x < -2.09 Then
        domflg = 15#
        t = Sqr(-x)
        zeta = -(2# * x * t / 3#)
        t = Sqr(t)
        k = sqpii / t
        z = 1# / zeta
        zz = z * z
        AFN = -0.131696323418332
        AFN = AFN * zz - 0.626456544431912
        AFN = AFN * zz - 0.693158036036933
        AFN = AFN * zz - 0.279779981545119
        AFN = AFN * zz - 0.04919001326095
        AFN = AFN * zz - 4.06265923594885E-03
        AFN = AFN * zz - 1.59276496239262E-04
        AFN = AFN * zz - 2.77649108155233E-06
        AFN = AFN * zz - 1.67787698489115E-08
        AFD = 1#
        AFD = AFD * zz + 13.3560420706553
        AFD = AFD * zz + 32.6825032795225
        AFD = AFD * zz + 26.73670409415
        AFD = AFD * zz + 9.1870740290726
        AFD = AFD * zz + 1.47529146771666
        AFD = AFD * zz + 0.115687173795188
        AFD = AFD * zz + 4.40291641615211E-03
        AFD = AFD * zz + 7.54720348287414E-05
        AFD = AFD * zz + 4.5185009297058E-07
        uf = 1# + zz * AFN / AFD
        AGN = 1.97339932091686E-02
        AGN = AGN * zz + 0.391103029615688
        AGN = AGN * zz + 1.06579897599596
        AGN = AGN * zz + 0.93916922981665
        AGN = AGN * zz + 0.351465656105548
        AGN = AGN * zz + 6.33888919628925E-02
        AGN = AGN * zz + 5.85804113048388E-03
        AGN = AGN * zz + 2.82851600836737E-04
        AGN = AGN * zz + 6.98793669997261E-06
        AGN = AGN * zz + 8.11789239554389E-08
        AGN = AGN * zz + 3.41551784765924E-10
        AGD = 1#
        AGD = AGD * zz + 9.30892908077442
        AGD = AGD * zz + 19.8352928718312
        AGD = AGD * zz + 15.5646628932865
        AGD = AGD * zz + 5.47686069422975
        AGD = AGD * zz + 0.954293611618962
        AGD = AGD * zz + 8.64580826352392E-02
        AGD = AGD * zz + 4.12656523824223E-03
        AGD = AGD * zz + 1.01259085116509E-04
        AGD = AGD * zz + 1.17166733214414E-06
        AGD = AGD * zz + 4.9183457006293E-09
        ug = z * AGN / AGD
        theta = zeta + 0.25 * PI()
        f = Sin(theta)
        g = Cos(theta)
        Ai = k * (f * uf - g * ug)
        Bi = k * (g * uf + f * ug)
        APFN = 0.185365624022536
        APFN = APFN * zz + 0.886712188052584
        APFN = APFN * zz + 0.987391981747399
        APFN = APFN * zz + 0.401241082318004
        APFN = APFN * zz + 7.10304926289631E-02
        APFN = APFN * zz + 5.90618657995662E-03
        APFN = APFN * zz + 2.33051409401777E-04
        APFN = APFN * zz + 4.08718778289035E-06
        APFN = APFN * zz + 2.48379932900442E-08
        APFD = 1#
        APFD = APFD * zz + 14.7345854687503
        APFD = APFD * zz + 37.542393343549
        APFD = APFD * zz + 31.4657751203046
        APFD = APFD * zz + 10.9969125207299
        APFD = APFD * zz + 1.78885054766999
        APFD = APFD * zz + 0.141733275753663
        APFD = APFD * zz + 5.44066067017226E-03
        APFD = APFD * zz + 9.39421290654511E-05
        APFD = APFD * zz + 5.65978713036027E-07
        uf = 1# + zz * APFN / APFD
        APGN = -3.55615429033082E-02
        APGN = APGN * zz - 0.637311518129436
        APGN = APGN * zz - 1.70856738884312
        APGN = APGN * zz - 1.50221872117317
        APGN = APGN * zz - 0.563606665822103
        APGN = APGN * zz - 0.102101031120217
        APGN = APGN * zz - 9.48396695961445E-03
        APGN = APGN * zz - 4.60325307486781E-04
        APGN = APGN * zz - 1.14300836484517E-05
        APGN = APGN * zz - 1.33415518685547E-07
        APGN = APGN * zz - 5.63803833958894E-10
        APGD = 1#
        APGD = APGD * zz + 9.8586580169613
        APGD = APGD * zz + 21.6401867356586
        APGD = APGD * zz + 17.3130776389749
        APGD = APGD * zz + 6.17872175280829
        APGD = APGD * zz + 1.08848694396321
        APGD = APGD * zz + 9.95005543440888E-02
        APGD = APGD * zz + 4.78468199683887E-03
        APGD = APGD * zz + 1.18159633322839E-04
        APGD = APGD * zz + 1.37480673554219E-06
        APGD = APGD * zz + 5.79912514929148E-09
        ug = z * APGN / APGD
        k = sqpii * t
        Aip = -(k * (g * uf + f * ug))
        Bip = k * (f * uf - g * ug)
        Exit Sub
    End If
    If x >= 2.09 Then
        domflg = 5#
        t = Sqr(x)
        zeta = 2# * x * t / 3#
        g = Exp(zeta)
        t = Sqr(t)
        k = 2# * t * g
        z = 1# / zeta
        AN = 0.346538101525629
        AN = AN * z + 12.0075952739646
        AN = AN * z + 76.2796053615235
        AN = AN * z + 168.089224934631
        AN = AN * z + 159.756391350164
        AN = AN * z + 70.5360906840444
        AN = AN * z + 14.026469116339
        AN = AN * z + 1#
        AD = 0.56759453263877
        AD = AD * z + 14.7562562584847
        AD = AD * z + 84.5138970141475
        AD = AD * z + 177.3180881454
        AD = AD * z + 164.23469287153
        AD = AD * z + 71.4778400825576
        AD = AD * z + 14.0959135607834
        AD = AD * z + 1#
        f = AN / AD
        Ai = sqpii * f / k
        k = -(0.5 * sqpii * t / g)
        APN = 0.613759184814036
        APN = APN * z + 14.7454670787755
        APN = APN * z + 82.0584123476061
        APN = APN * z + 171.184781360976
        APN = APN * z + 159.317847137142
        APN = APN * z + 69.9778599330103
        APN = APN * z + 13.9470856980482
        APN = APN * z + 1#
        APD = 0.334203677749737
        APD = APD * z + 11.1810297306158
        APD = APD * z + 71.172735214786
        APD = APD * z + 158.778084372838
        APD = APD * z + 153.206427475809
        APD = APD * z + 68.675230459278
        APD = APD * z + 13.8498634758259
        APD = APD * z + 1#
        f = APN / APD
        Aip = f * k
        If x > 8.3203353 Then
            BN16 = -0.253240795869364
            BN16 = BN16 * z + 0.575285167332467
            BN16 = BN16 * z - 0.329907036873225
            BN16 = BN16 * z + 0.06444040689482
            BN16 = BN16 * z - 3.82519546641337E-03
            BD16 = 1#
            BD16 = BD16 * z - 7.15685095054035
            BD16 = BD16 * z + 10.6039580715665
            BD16 = BD16 * z - 5.23246636471251
            BD16 = BD16 * z + 0.957395864378384
            BD16 = BD16 * z - 0.055082814716355
            f = z * BN16 / BD16
            k = sqpii * g
            Bi = k * (1# + f) / t
            BPPN = 0.465461162774652
            BPPN = BPPN * z - 1.08992173800494
            BPPN = BPPN * z + 0.638800117371828
            BPPN = BPPN * z - 0.126844349553103
            BPPN = BPPN * z + 7.6248784434211E-03
            BPPD = 1#
            BPPD = BPPD * z - 8.70622787633159
            BPPD = BPPD * z + 13.8993162704553
            BPPD = BPPD * z - 7.14116144616431
            BPPD = BPPD * z + 1.34008595960681
            BPPD = BPPD * z - 7.84273211323342E-02
            f = z * BPPN / BPPD
            Bip = k * t * (1# + f)
            Exit Sub
        End If
    End If
    f = 1#
    g = x
    t = 1#
    uf = 1#
    ug = x
    k = 1#
    z = x * x * x
    Do While t > MachineEpsilon
        uf = uf * z
        k = k + 1#
        uf = uf / k
        ug = ug * z
        k = k + 1#
        ug = ug / k
        uf = uf / k
        f = f + uf
        k = k + 1#
        ug = ug / k
        g = g + ug
        t = Abs(uf / f)
    Loop
    uf = c1 * f
    ug = c2 * g
    If domflg Mod 2# = 0# Then
        Ai = uf - ug
    End If
    If domflg \ 2# Mod 2# = 0# Then
        Bi = sqrt3 * (uf + ug)
    End If
    k = 4#
    uf = x * x / 2#
    ug = z / 3#
    f = uf
    g = 1# + ug
    uf = uf / 3#
    t = 1#
    Do While t > MachineEpsilon
        uf = uf * z
        ug = ug / k
        k = k + 1#
        ug = ug * z
        uf = uf / k
        f = f + uf
        k = k + 1#
        ug = ug / k
        uf = uf / k
        g = g + ug
        k = k + 1#
        t = Abs(ug / g)
    Loop
    uf = c1 * f
    ug = c2 * g
    If domflg \ 4# Mod 2# = 0# Then
        Aip = uf - ug
    End If
    If domflg \ 8# Mod 2# = 0# Then
        Bip = sqrt3 * (uf + ug)
    End If
End Sub