Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

made topolar() more accurate, added fuzztest for Calc_R() #29

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
77 changes: 53 additions & 24 deletions src/ComplexHuff/Complex.huff
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,7 @@
#define macro TO_POLAR() = takes(2) returns(2) {
// INPUT STACK => [RE(a),IM(a)]


dup2 //[i,r,i]
dup2 //[r,i,r,i]
CALC_R() //[r,ra,ia]
Expand All @@ -116,28 +117,28 @@
dup1
[X3]
eq case2 jumpi
[X3]
swap1
sub //[x1-1e18]
0x02
swap1
sdiv // [x1_new]
dup1
0x02
mul
[X13]
swap1
sdiv
0x02
swap1
exp
0x04
swap1
sdiv
0x01
0x00
sub
mul //[x2_new,x1_new]
[X3] //[1e18,x1,r]
swap1 //[x1,1e18,r]
sub //[x1-1e18,r]
0x02 //[2,x1-1e18,r]
swap1 //[x1-1e18,2,r]
sdiv //[x1_new=x1-1e18/2,r]
dup1 //[x1_new,x1_new,r]
0x02 //[2,x1_new,x1_new,r]
mul //[2*x1_new,x1_new,r]
[X13] //[1e9,...]
swap1 //[2*x1_new,1e9,..]
sdiv //[2*x1_new/1e9,..]
0x02 //[2,...]
swap1 //[a=2*x1_new/1e9,2,..]
exp //[a**2,..]
0x04 //[4,...]
swap1 //[a**2,4,..]
sdiv //[a**2/4,..]
0x01 //[1,...]
0x00 //[0,1,...]
sub //[-1,...]
mul //[x2_new=]
dup2
0x02
mul
Expand Down Expand Up @@ -167,7 +168,7 @@
0x03 // [3,a**3,x1,r]
swap1 // [a**3,3,..]
sdiv // [a**3/3 , x, r]
0x01 // [1,a**3/3,..]
0x01 // [1,a**3/3,..]
0x00 // [0,1,a**3/3,..]
sub // [-1,a**3/3,..]
mul // [x2=-(a**3)/3,x1,r]
Expand All @@ -180,8 +181,36 @@
exp // [b**5,..]
0x05
swap1
sdiv // [x3,x2,x1,r]
sdiv // [x3,x2,x1,r]
add // [x3+x2,x1,r]
dup2 // [x1,x3+x2,x1,r]
[X16] // [X16,x1,..]
swap1 // [x1,X16,...]
sdiv // [c=x1/X16,x3+x2,x1,r]
0x07 // [7,c,...]
swap1 // [c,7,...]
exp // [c**7,...]
0x07 // [7,c**7,...]
swap1 // [c**7,7,...]
sdiv // [c**7/7,...]
0x01 // [1,....]
0x00 // [0,1,....]
sub // [-1,....]
mul // [-c**7/7,...]
add // [x4+x3+x2,x1,r]
dup2 // [x1,...]
[X17] // [X17,x1,...]
swap1 // [x1,X17,...]
sdiv // [x1/X17,...]
0x08 // [8,x1/X17,...]
swap1 // [x1/X17,8,...]
exp // [d**8,...]
0x960 // [2400,d**8,...]
mul // [2400*d**8,...]
0x9d80 // [40320,2400*d**8,...]
swap1 // [2400*d**8,40320,...]
sdiv // [x5,x4+x3+x2,x1,r] // x5 is calculated by error approximation and hit and trial to make the graph similar to arctanx
add
add // [theta,r]
swap1
finish jump
Expand Down
7 changes: 6 additions & 1 deletion src/ComplexHuff/Constants.huff
Original file line number Diff line number Diff line change
Expand Up @@ -16,4 +16,9 @@
#define constant HALF_SCALE = 0x6F05B59D3B20000 //5e17
#define constant X13 = 0x3B9ACA00 // 1e9
#define constant X14 = 0x1CC2C05DBC54 // 1e(13.5)
#define constant X15 = 0x2386F26FC10000 // 1e16
#define constant X15 = 0x2386F26FC10000 // 1e16
#define constant X16 = 0x987e5c9eb0764 // 1e(18*6/7)
#define constant X17 = 0x13fa76ed4e1003 // 1e(18*7/8)



63 changes: 33 additions & 30 deletions test/Complex.t.sol
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ contract ComplexTest is Test {
function testToPolar() public {
(int256 r, int256 t) = complex.toPolar(3, 4);
assertEq(r, 5);
assertEq((t * 100) / scale, 65);
assertEq(t / 1e15, 643); // arctan(0.75)=0.643501108793 approximated to 3 decimal places
}

function testFromPolar() public {
Expand Down Expand Up @@ -78,7 +78,7 @@ contract ComplexTest is Test {
function testLnZ() public {
(int256 r, int256 i) = complex.ln(30 * scale, 40 * scale);
assertEq((r * 100) / scale, 391); // ln(50) = 3.912..
assertEq((i * 100) / scale, 65);
assertEq((i * 100) / scale, 64);
}

function testAtan2() public {
Expand Down Expand Up @@ -254,34 +254,37 @@ contract ComplexTest is Test {
assertEq(Ri, Rhi);
}

// function testCalc_RFuzz(int256 ar, int256 ai, int256 br, int256 bi, int256 k) public {
// ai = bound(ai, -1e9, 1e9);
// bi = bound(bi, -1e9, 1e9);
// ar = bound(ar, -1e9, 1e9);
// br = bound(br, -1e9, 1e9);
// k = bound(k, -1e9, 1e9);

// (int256 mag1r,) = (complex.mulz(-ai, ai, ar, ar)); // ar^2 + ai^2
// int256 mag1 = PRBMathSD59x18.sqrt(mag1r); // (ar^2+ai^2)^0.5
// uint256 R1 = (complex.calcR(ai, ar)); // magnitude(A)
// uint256 R2 = (complex.calcR(bi, br)); // magnitude(B)
// (int256 A_Br, int256 A_Bi) = (complex.addz(bi, ai, br, ai)); // A+B
// uint256 R3 = (complex.calcR(A_Bi, A_Br)); // magnitude(A+B)
// uint256 R4 = (complex.calcR(k * ai, k * ar)); // magnitude(k*A)
// uint256 magk = (complex.calcR(0, k));
// uint256 _R1 = (complex.calcR(-ai, ar));

// // Test by comparing
// assertEq(uint256(mag1) / 1e9, R1);

// // Test by property
// // mag(A+B)<=mag(A)+mag(B)
// assert(R3 <= (R1 + R2));
// // mag(kA)=kmag(A)
// assertEq(R4,(magk)*R1);
// // mag(A)=mag(A')
// assertEq(R1,_R1);
// }
function testCalc_RFuzz(int256 ar, int256 ai, int256 br, int256 bi, int256 k) public {
ai = bound(ai, -1e9, 1e9);
bi = bound(bi, -1e9, 1e9);
ar = bound(ar, -1e9, 1e9);
br = bound(br, -1e9, 1e9);
k = bound(k, -1e9, 1e9);
vm.assume(ai != 0);
vm.assume(bi != 0);
vm.assume(ar != 0);
vm.assume(br != 0);

(int256 mag1r,) = (complex.mulz(-ai * scale, ai * scale, ar * scale, ar * scale)); // ar^2 + ai^2
mag1r = PRBMathSD59x18.sqrt(mag1r); // (ar^2+ai^2)^0.5
uint256 R1 = (complex.calcR(ai * scale, ar * scale)); // magnitude(A)
// Test by comparing
assertEq(uint256(mag1r) / 1e9, R1);

// uint256 R2 = (complex.calcR(bi*scale, br*scale)); // magnitude(B)
// (int256 A_Br, int256 A_Bi) = (complex.addz(bi*scale, ai*scale, br*scale, ai*scale)); // A+B
// uint256 R3 = (complex.calcR(A_Bi, A_Br)); // magnitude(A+B)
// uint256 R4 = (complex.calcR(k * ai*scale, k * ar*scale)); // magnitude(k*A)
// uint256 magk = (complex.calcR(0, k));
// uint256 _R1 = (complex.calcR(-ai*scale, ar*scale));
// // Test by property
// // mag(A+B)<=mag(A)+mag(B)
// assertGe(R1 + R2 ,R3 ); //This test is failing
// // // mag(kA)=mag(k)*mag(A)
// // assertEq(int(R4)/scale,int((magk)*R1)/scale);
// // // mag(A)=mag(A')
// // assertEq(R1,_R1);
}
}

interface WRAPPER {
Expand Down
Loading