Files
ground-programs/packages/math/math.grnd

583 lines
11 KiB
Plaintext

fun -double !mod -double &in -double &modulus
lesser $in $modulus &store
not $store &store
if $store %downloop
lesser $in 0 &store
if $store %uploop
jump %end
@downloop
subtract $in $modulus &in
lesser $in $modulus &store
not $store &store
if $store %downloop
jump %end
@uploop
add $in $modulus &in
lesser $in 0 &store
if $store %uploop
jump %end
@end
return $in
endfun
fun -double !pi
return 3.141592653589793
endfun
fun -double !e
return 2.718281828459045
endfun
fun -double !intexp -double &base -int &power
lesser $power 1 &store
set &ans 1
if $store %end
@loop
subtract $power 1 &power
multiply $ans $base &ans
lesser $power 1 &store
if $store %end
jump %loop
@end
return $ans
endfun
fun -int !factorial -int &in
set &idx $in
set &ans 1
lesser $idx 0 &store
if $store %error
@loop
equal $idx 0 &store
if $store %end
multiply $ans $idx &ans
subtract $idx 1 &idx
jump %loop
@end
return $ans
@error
error "Cannot find factorial of a negative integer"
endfun
fun -double !cos -double &in
# Get input modulo 2π
call !math:pi &pi
multiply $pi 2 &tau
pusharg $in $tau
call !math:mod &in
set &idx 0
set &ans 0
# Check if $in < π
greater $in $pi &store
if $store %cos2
jump %cos1
@cos1
# Taylor series at x=π/2
multiply $idx 2 &2k+1
add $2k+1 1 &2k+1
# Get (x-π/2)^(2k+1)
divide $pi 2 &store
subtract $in $store &store
pusharg $store
tostring $2k+1 &2k+1
stoi $2k+1 &2k+1
pusharg $2k+1
call !math:intexp &exp
# Get (2k+1)!
pusharg $2k+1
call !math:factorial &fact
# Divide
divide $exp $fact &tempans
# Add result to $ans
tostring $idx &idx
stod $idx &idx
pusharg $idx
stod "2" &store
pusharg $store
call !math:mod &store
equal $store 0 &store
if $store %subtract
jump %add
@add
add $ans $tempans &ans
jump %increment
@subtract
subtract $ans $tempans &ans
jump %increment
@increment
add $idx 1 &idx
equal $idx 5 &store
if $store %end
jump %cos1
@cos2
# cos(x)=cos(2π-x)
subtract $tau $in &in
jump %cos1
@end
return $ans
endfun
fun -double !sin -double &in
call !math:pi &halfpi
divide $halfpi 2 &halfpi
subtract $halfpi $in &in
pusharg $in
call !math:cos &store
return $store
endfun
fun -double !tan -double &in
pusharg $in
call !math:sin &store
pusharg $in
call !math:cos &store2
divide $store $store2 &store
return $store
endfun
fun -double !asin -double &in
pusharg $in
call !math:abs &store
greater $store 1 &store
if $store %error
# asin(a) is the root of sin(x)-a=0, -pi/2 <= x <= pi/2
# Using Newton's Method: x_k+1 = x_k - (sin(x_k)-a)/cos(x_k)
stod "0" &ans
@loop
pusharg $ans
call !math:sin &sin
subtract $sin $in &sin
pusharg $ans
call !math:cos &cos
divide $sin $cos &store
subtract $ans $store &ans
equal $store 0 &store
if $store %return
jump %loop
@return
return $ans
@error
error "Cannot find sin of a double with magnetude greater than 1" "mathError"
endfun
fun -double !acos -double &in
# asin(x) + acos(x) = pi/2
pusharg $in
call !math:asin &store
call !math:pi &pi
divide $pi 2 &pi
subtract $pi $store &store
return $store
endfun
fun -double !atan -double &in
# atan(a) is the root of tan(x)-a=0, -pi/2 <= x <= pi/2
# Using Newton's Method: x_k+1 = x_k - (tan(x_k)-a)/sec^2(x_k) = x_k - (tan(x_k)-a)cos^2(x_k)
stod "0" &ans
@loop
pusharg $ans
call !math:tan &num
subtract $num $in &num
pusharg $ans
call !math:cos &dom
multiply $dom $dom &dom
multiply $num $dom &store
subtract $ans $store &ans
pusharg $store
call !math:abs &store
lesser $store 0.000000000001 &store
if $store %return
jump %loop
@return
# Make answer in range [-pi/2, pi/2]
call !math:pi &pi
pusharg $ans
pusharg $pi
call !math:mod &ans
divide $pi 2 &pi
greater $ans $pi &store
if $store %negans
return $ans
@negans
call !math:pi &pi
subtract $ans $pi &ans
return $ans
endfun
fun -double !atan2 -double &y -double &x
call !math:pi &pi
divide $pi 2 &pi/2
# catch "mathError" &divideByZero
set &divideByZero true
divide $y $x &ratio
if $divideByZero %continue
jump %error
@continue
# When x > 0, return atan(y/x). When x < 0, return pi+atan(y/x) mod 2pi
greater $x 0 &store
if $store %atan
# When x < 0, check sign of y
lesser $y 0 &store
if $store %quadrant3
pusharg $ratio
call !math:atan &store
add $store $pi &store
return $store
@quadrant3
pusharg $ratio
call !math:atan &store
subtract $store $pi &store
return $store
@atan
pusharg $ratio
call !math:atan &store
return $store
@error
greater $y 0 &store
if $store %90deg
lesser $y 0 &store
if $store %neg90deg
error "Cannot find atan of (0,0)" "mathError"
divide $pi 2 &pi/2
@90deg
return $pi/2
@neg90deg
multiply $pi/2 -1 &store
return $store
endfun
fun -double !random -int &seed
set &m 2147483648
set &a 1103515245
set &c 12345
multiply $seed $a &store
add $store $c &store
tostring $store &store
stod $store &store
pusharg $store
pusharg $m
call !math:mod &store
divide $store 2147483648 &store
return $store
endfun
fun -double !itod -int &int
tostring $int &int
stod $int &int
return $int
endfun
fun -int !dtoi -double &double
tostring $double &double
stoi $double &double
return $double
endfun
fun -double !sqrt -double &in
# Finding sqrt(a) is the same as finding the roots of
# f(x) = x^2 - a
# Using Newton's method: x_(k+1) = x_k - f(x_k)/f'(x_k) = x_k - ((x_k)^2-a)/2(x_k)
lesser $in 0 &store
if $store %error
set &ans 2
@loop
multiply $ans $ans &num
subtract $num $in &num
multiply $ans 2 &dom
divide $num $dom &store
equal $store 0 &store2
if $store2 %end
subtract $ans $store &ans
jump %loop
@end
return $ans
@error
error "Cannot find square root of a negative double" "mathError"
endfun
fun -double !abs -double &in
lesser $in 0 &store
if $store %negative
return $in
@negative
multiply $in -1 &in
return $in
endfun
fun -double !exp -double &in
call !math:e &e
lesser $in 0 &store
if $store %negative
# Find $in rounded down
pusharg $in
stod "1.00" &store
pusharg $store
call !math:mod &floorin
subtract $in $floorin &floorin
pusharg $floorin
call !math:dtoi &floorin
# e^(a+x)=e^a*e^x
pusharg $e
pusharg $floorin
call !math:intexp &e^a
# Find e^($in - $floorin)
subtract $in $floorin &x
# e^x = 1 + x + x^2/2 + x^3/6 + ... + x^k/k!
set &idx 0
set &e^x 0
@loop
equal $idx 10 &store
if $store %return
pusharg $x
pusharg $idx
call !math:intexp &store
pusharg $idx
call !math:factorial &fact
divide $store $fact &store
add $e^x $store &e^x
add $idx 1 &idx
jump %loop
@return
multiply $e^a $e^x &ans
return $ans
@negative
# e^-x = 1/e^x
multiply $in -1 &in
pusharg $in
call !math:exp &store
divide 1 $store &store
return $store
endfun
fun -double !ln -double &in
greater $in 0 &store
if $store %calculate
jump %error
@calculate
equal $in 1 &store
if $store %zero
# We will use Newton's method.
# ln(a) is the root of e^x-a
# Therefore x_k+1 = x_k - (e^x_k - a)/(e^x_k)
set &ans 0.9
@loop
pusharg $ans
call !math:exp &e^x
subtract $e^x $in &store
divide $store $e^x &store
subtract $ans $store &ans
pusharg $store
call !math:abs &store
lesser $store 0.000000000001 &store
if $store %return
jump %loop
@return
return $ans
@zero
stod "0" &zero
return $zero
@error
error "Cannot find ln of a non-positive double" "mathError"
endfun
fun -double !rpn -list &rpn
use "lists"
# List format: [3, 2, 4, '+', '-'] returns -3
setlist *stack
getlistsize *rpn &len
set &idx 0
@loop
equal $idx $len &store
if $store %return
# If type is a number, add it to the stack
getlistat *rpn $idx &store
gettype $store &store
equal $store "double" &bool
if $bool %pushstack
equal $store "int" &bool
if $bool %pushstack
getlistat *rpn $idx &store
equal $store '+' &bool
if $bool %add
equal $store '-' &bool
if $bool %subtract
equal $store '*' &bool
if $bool %multiply
equal $store '/' &bool
if $bool %divide
jump %increment
@add
getlistsize *stack &store
subtract $store 2 &store
getlistat *stack $store &list1
add $store 1 &store
getlistat *stack $store &list2
add $list1 $list2 &list1
subtract $store 1 &store
setlistat *stack $store $list1
pusharg *stack
call !lists:listremove &stack
jump %increment
@subtract
getlistsize *stack &store
subtract $store 2 &store
getlistat *stack $store &list1
add $store 1 &store
getlistat *stack $store &list2
subtract $list1 $list2 &list1
subtract $store 1 &store
setlistat *stack $store $list1
pusharg *stack
call !lists:listremove &stack
jump %increment
@multiply
getlistsize *stack &store
subtract $store 2 &store
getlistat *stack $store &list1
add $store 1 &store
getlistat *stack $store &list2
multiply $list1 $list2 &list1
subtract $store 1 &store
setlistat *stack $store $list1
pusharg *stack
call !lists:listremove &stack
jump %increment
@divide
getlistsize *stack &store
subtract $store 2 &store
getlistat *stack $store &list1
add $store 1 &store
getlistat *stack $store &list2
divide $list1 $list2 &list1
subtract $store 1 &store
setlistat *stack $store $list1
pusharg *stack
call !lists:listremove &stack
jump %increment
@pushstack
getlistat *rpn $idx &store
listappend *stack $store
jump %increment
@increment
add $idx 1 &idx
jump %loop
@return
getlistsize *stack &store
greater $store 1 &store
if $store %error
getlistat *stack 0 &store
return $store
@error
error "Stack for !fpn has multiple contents" "rangeError"
endfun