/// @title Trigonometry
/// @notice SPDX-License-Identifier: MIT
/// @author mds1 <https://github.com/mds1>
/// @author clabby <https://github.com/clabby>
/// @notice Basic trigonometry functions where inputs and outputs are integers.
/// Inputs are specified in radians scaled by 1e18, and similarly outputs are scaled by 1e18.
/// @notice Adapted from mds1 (https://github.com/mds1/solidity-trigonometry/blob/main/src/Trigonometry.sol)
////////////////////////////////////////////////////////////////
// CONSTANTS //
////////////////////////////////////////////////////////////////
// Table index into the trigonometric table
#define constant INDEX_WIDTH = 0x08
// Interpolation between successive entries in the table
#define constant INTERP_WIDTH = 0x10
#define constant INDEX_OFFSET = 0x14
#define constant INTERP_OFFSET = 0x04
#define constant ANGLES_IN_CYCLE = 0x40000000
#define constant QUADRANT_HIGH_MASK = 0x20000000
#define constant QUADRANT_LOW_MASK = 0x10000000
#define constant SINE_TABLE_SIZE = 0x100
// Pi as an 18 decimal value, which is plenty of accuracy: "For JPL's highest accuracy calculations, which are for
// interplanetary navigation, we use 3.141592653589793: https://www.jpl.nasa.gov/edu/news/2016/3/16/how-many-decimals-of-pi-do-we-really-need/
#define constant PI = 0x2B992DDFA23249D6
#define constant TWO_PI = 0x57325BBF446493AC
#define constant PI_OVER_TWO = 0x15CC96EFD1192500
// The constant sine lookup table was generated by generate_trigonometry.py. We must use a constant
// bytes array because constant arrays are not supported in Solidity. Each entry in the lookup
// table is 4 bytes. Since we're using 32-bit parameters for the lookup table, we get a table size
// of 2^(32/4) + 1 = 257, where the first and last entries are equivalent (hence the table size of
// 256 defined above)
#define constant ENTRY_BYTES = 0x04 // each entry in the lookup table is 4 bytes
#define constant ENTRY_MASK = 0xFFFFFFFF // mask used to cast bytes32 -> lookup table entry
#define table SIN_LUT {
0x00000000000000000000000000000000000000000000000000000000000001000000000000c90f8801921d20025b26d703242abf03ed26e604b6195d057f00350647d97c0710a34507d95b9e08a2009a096a90490a3308bc0afb68050bc3ac350c8bd35e0d53db920e1bc2e40ee387660fab272b1072a0481139f0cf120116d512c8106e138edbb1145576b1151bdf8515e2144416a81305176dd9de183366e818f8b83c19bdcbf31a82a0251b4732ef1c0b826a1ccf8cb31d934fe51e56ca1e1f19f97b1fdcdc1b209f701c2161b39f2223a4c522e541af23a6887e2467775725280c5d25e845b626a8218527679df42826b92828e5714a29a3c4852a61b1012b1f34eb2bdc4e6f2c98fbba2d553afb2e110a622ecc681e2f8752623041c76030fbc54d31b54a5d326e54c73326e2c233def2873496824f354d905636041ad936ba2013376f9e46382493b038d8fe93398cdd323a402dd13af2eeb73ba51e293c56ba703d07c1d53db832a53e680b2c3f1749b73fc5ec974073f21d4121589a41ce1e64427a41d04325c13543d09aec447acd50452456bc45cd358f46756827471cece647c3c22e4869e664490f57ee49b415334a581c9d4afb6c974b9e038f4c3fdff34ce100344d8162c34e2106174ebfe8a44f5e08e24ffb654c5097fc5e5133cc9451ced46e5269126e53028517539b2aef5433027d54ca0a4a556040e255f5a4d2568a34a9571deef957b0d2555842dd5458d40e8c5964649759f3de125a8279995b1035ce5b9d11535c290acc5cb420df5d3e52365dc79d7b5e50015d5ed77c895f5e0db25fe3b38d60686cce60ec382f616f146b61f1003e6271fa6862f201ac637114cc63ef328f646c59bf64e889256563bf9165ddfbd266573cbb66cf811f6746c7d767bd0fbc683257aa68a69e806919e31f698c246b69fd614a6a6d98a36adcc9646b4af2786bb812d06c24295f6c8f351b6cf934fb6d6227f96dca0d146e30e3496e96a99c6efb5f116f5f02b16fc1938470231099708378fe70e2cbc571410804719e2cd171fa394872552c8472af05a67307c3cf735f662573b5ebd0740b53fa745f9dd074b2c8837504d3447555bd4b75a585ce75f42c0a7641af3c768e0ea576d9498877235f2c776c4eda77b417df77fab988784033287884841378c7aba17909a92c794a7c11798a23b079c89f6d7a05eeac7a4210d87a7d055a7ab6cba37aef63237b26cb4e7b5d039d7b920b887bc5e28f7bf8882f7c29fbed7c5a3d4f7c894bdd7cb727237ce3ceb17d0f42177d3980eb7d628ac57d8a5f3f7db0fdf77dd6668e7dfa98a77e1d93e97e3f57fe7e5fe4927e7f39567e9d55fb7eba3a387ed5e5c57ef0585f7f0991c37f2191b37f3857f57f4de4507f62368e7f754e7f7f872bf27f97cebc7fa736b37fb563b27fc255957fce0c3d7fd8878d7fe1c76a7fe9cbbf7ff094777ff621817ffa72d07ffd88597fff62157fffffff
}
#define constant WAD = 0x0de0b6b3a7640000
#define constant I32_MAX = 0x7fffffff
////////////////////////////////////////////////////////////////
// BASIC TRIG //
////////////////////////////////////////////////////////////////
/// @notice Return the sine of a value, specified in radians scaled by 1e18
/// @dev This algorithm for converting sine only uses integer values, and it works by dividing the
/// circle into 30 bit angles, i.e. there are 1,073,741,824 (2^30) angle units, instead of the
/// standard 360 degrees (2pi radians). From there, we get an output in range -2,147,483,647 to
/// 2,147,483,647, (which is the max value of an int32) which is then converted back to the standard
/// range of -1 to 1, again scaled by 1e18
/// @param _angle Angle to convert
/// @return Result scaled by 1e18
#define macro SIN() = takes (1) returns (1) {
// Input stack: [angle]
// Convert angle from from arbitrary radian value (range of 0 to 2pi) to the algorithm's range
// of 0 to 1,073,741,824
[TWO_PI] dup1 // [TWO_PI, TWO_PI, angle]
swap2 // [angle, TWO_PI, TWO_PI]
mod // [angle % TWO_PI, TWO_PI]
[ANGLES_IN_CYCLE] mul // [ANGLES_IN_CYCLE * (angle % TWO_PI), TWO_PI]
div // [ANGLES_IN_CYCLE * (angle % TWO_PI) / TWO_PI]
// [angle]
// Apply a mask on an integer to extract a certain number of bits, where angle is the integer
// whose bits we want to get, the width is the width of the bits (in bits) we want to extract,
// and the offset is the offset of the bits (in bits) we want to extract. The result is an
// integer containing _width bits of _value starting at the offset bit
dup1 // [angle, angle]
[INTERP_OFFSET] shr // [angle >> INTERP_OFFSET, angle]
0x01 dup1 // [0x01, 0x01, angle >> INTERP_OFFSET, angle]
[INTERP_WIDTH] shl // [0x01 << INTERP_WIDTH, 0x01, angle >> INTERP_OFFSET, angle]
sub // [(0x01 << INTERP_WIDTH) - 0x01, angle >> INTERP_OFFSET, angle]
and // [interp, angle]
dup2 // [angle, interp, angle]
[INDEX_OFFSET] shr // [angle >> INDEX_OFFSET, interp, angle]
0x01 dup1 // [0x01, 0x01, angle >> INDEX_OFFSET, interp, angle]
[INDEX_WIDTH] shl // [0x01 << INDEX_WIDTH, 0x01, angle >> INDEX_OFFSET, interp, angle]
sub // [(0x01 << INDEX_WIDTH) - 0x01, angle >> INDEX_OFFSET, interp, angle]
and // [index, interp, angle]
// The lookup table only contains data for one quadrant (since sin is symmetric around both
// axes), so here we figure out which quadrant we're in, then we lookup the values in the
// table then modify values accordingly
dup3 // [angle, index, interp, angle]
[QUADRANT_LOW_MASK] // [QUADRANT_LOW_MASK, angle, index, interp, angle]
and // [QUADRANT_LOW_MASK & angle, index, interp, angle]
iszero // [is_odd_quadrant, index, interp, angle]
dup4 // [angle, is_odd_quadrant, index, interp, angle]
[QUADRANT_HIGH_MASK] // [QUADRANT_HIGH_MASK, angle, is_odd_quadrant, index, interp, angle]
and // [QUADRANT_HIGH_MASK & angle, is_odd_quadrant, index, interp, angle]
iszero iszero // [is_negative_quadrant, is_odd_quadrant, index, interp, angle]
// Jump past updating the index if `is_odd_quadrant` is true
dup2 is_odd_q jumpi // [is_negative_quadrant, is_odd_quadrant, index, interp, angle]
dup3 // [index, is_negative_quadrant, is_odd_quadrant, index, interp, angle]
0x01 // [0x01, index, is_negative_quadrant, is_odd_quadrant, index, interp, angle]
[SINE_TABLE_SIZE] // [SINE_TABLE_SIZE, 0x01, index, is_negative_quadrant, is_odd_quadrant, index, interp, angle]
sub sub // [SINE_TABLE_SIZE - 0x01 - index, is_negative_quadrant, is_odd_quadrant, index, interp, angle]
swap3 pop // [is_negative_quadrant, is_odd_quadrant, index, interp, angle]
is_odd_q:
// We are looking for two consecutive indices in our lookup table
// Since EVM is left aligned, to read n bytes of data from idx i, we must read from `i * data_len` + `n`
// therefore, to read two entries of size entry_bytes `index * entry_bytes` + `entry_bytes * 2`
swap2 // [index, is_odd_quadrant, is_negative_quadrant, interp, angle]
0x02 add // [index + 0x02, is_odd_quadrant, is_negative_quadrant, interp, angle]
[ENTRY_BYTES] mul // [offset1_2, is_odd_quadrant, is_negative_quadrant, interp, angle]
// This following snippet will function for any entry_bytes <= 15
__tablestart(SIN_LUT) // [sin_table_start, offset1_2, is_odd_quadrant, is_negative_quadrant, interp, angle]
add // [sin_table_start + offset1_2, is_odd_quadrant, is_negative_quadrant, interp, angle]
0x20 swap1 // [sin_table_start + offset1_2, 0x20, is_odd_quadrant, is_negative_quadrant, interp, angle]
0x00 codecopy // [is_odd_quadrant, is_negative_quadrant, interp, angle]
0x00 mload // [x1_2, is_odd_quadrant, is_negative_quadrant, interp, angle]
// We now read the last two numbers of size entry_bytes from x1_2
// in example: entry_bytes = 4; x1_2 = 0x00...12345678abcdefgh
// therefore: entry_mask = 0xFFFFFFFF
// 0x00...12345678abcdefgh >> 8*4 = 0x00...12345678
// 0x00...12345678 & 0xFFFFFFFF = 0x12345678
dup1 0x20 shr // [x1_2 >> 0x20, x1_2, is_odd_quadrant, is_negative_quadrant, interp, angle]
[ENTRY_MASK] and // [x1, x1_2, is_odd_quadrant, is_negative_quadrant, interp, angle]
// 0x00...12345678abcdefgh & 0xFFFFFFFF = 0xabcdefgh
swap1 [ENTRY_MASK] // [ENTRY_MASK, x1_2, x1, is_odd_quadrant, is_negative_quadrant, interp, angle]
and // [x2, x1, is_odd_quadrant, is_negative_quadrant, interp, angle]
// Approximate angle by interpolating in the table, accounting for the quadrant
dup2 dup2 sub // [x2 - x1, x2, x1, is_odd_quadrant, is_negative_quadrant, interp, angle]
dup6 mul // [interp * (x2 - x1), x2, x1, is_odd_quadrant, is_negative_quadrant, interp, angle]
[INTERP_WIDTH] shr // [(interp * (x2 - x1)) >> INTERP_WIDTH, x2, x1, is_odd_quadrant, is_negative_quadrant, interp, angle]
swap3 // [is_odd_quadrant, x2, x1, (interp * (x2 - x1)) >> INTERP_WIDTH, is_negative_quadrant, interp, angle]
is_odd jumpi // [x2, x1, (interp * (x2 - x1)) >> INTERP_WIDTH, is_negative_quadrant, interp, angle]
is_even:
swap1 pop // [x2, (interp * (x2 - x1)) >> INTERP_WIDTH, is_negative_quadrant, interp, angle]
sub // [x2 - (interp * (x2 - x1)) >> INTERP_WIDTH, is_negative_quadrant, interp, angle]
finish jump
is_odd:
pop add // [x1 + (interp * (x2 - x1)) >> INTERP_WIDTH, is_negative_quadrant, interp, angle]
finish:
// Received stack: [sine, is_negative_quadrant, interp, angle]
swap1 // [is_negative_quadrant, sine, interp, angle]
iszero cont jumpi // [sine, interp, angle]
// Negative 1
0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
mul // [-sine, interp, angle]
cont:
// Bring result from the range of -2,147,483,647 through 2,147,483,647 to -1e18 through 1e18.
// This can never overflow because sine is bounded by the above values
[WAD] mul // [1e18 * sine, interp, angle]
[I32_MAX] // [I32_MAX, 1e18 * sine, interp, angle]
swap1 sdiv // [1e18 * sine / I32_MAX, interp, angle]
// Clean stack:
swap2 pop pop // [1e18 * sine / I32_MAX]
// Return stack: [1e18 * sine / I32_MAX]
}
/// @notice Return the cosine of a value, specified in radians scaled by 1e18
/// @dev This is identical to the sin() method, and just computes the value by delegating to the
/// sin() method using the identity cos(x) = sin(x + pi/2)
/// @dev Overflow when `angle + PI_OVER_TWO > type(uint256).max` is ok, results are still accurate
/// @param _angle Angle to convert
/// @return Result scaled by 1e18
#define macro COS() = takes (1) returns (1) {
// Input stack: [angle]
[PI_OVER_TWO] add // [angle + pi/2]
SIN() // [result]
// Return stack: [result]
}