From 4dc2d4cfe511a798b30826a10bb0a121f0aabc39 Mon Sep 17 00:00:00 2001 From: Adrian Mariano Date: Fri, 12 Mar 2021 20:36:34 -0500 Subject: [PATCH] Fix bug in skin, add some complex stuff to math and tests --- arrays.scad | 13 +++++ math.scad | 114 ++++++++++++++++++++++++++++++++++++------- skin.scad | 2 +- tests/test_math.scad | 75 +++++++++++++++++++++++++--- 4 files changed, 177 insertions(+), 27 deletions(-) diff --git a/arrays.scad b/arrays.scad index cce3552..04e9a89 100644 --- a/arrays.scad +++ b/arrays.scad @@ -1682,4 +1682,17 @@ function transpose(arr, reverse=false) = arr; +// Function: is_matrix_symmetric() +// Usage: +// b = is_matrix_symmetric(A,) +// Description: +// Returns true if the input matrix is symmetric, meaning it equals its transpose. +// Matrix should have numerical entries. +// Arguments: +// A = matrix to test +// eps = epsilon for comparing equality. Default: 1e-12 +function is_matrix_symmetric(A,eps=1e-12) = + approx(A,transpose(A)); + + // vim: expandtab tabstop=4 shiftwidth=4 softtabstop=4 nowrap diff --git a/math.scad b/math.scad index e73c088..79d5ffa 100644 --- a/math.scad +++ b/math.scad @@ -1465,19 +1465,52 @@ function deriv3(data, h=1, closed=false) = // Section: Complex Numbers + +// Function: complex() +// Usage: +// z = complex(list) +// Description: +// Converts a real valued number, vector or matrix into its complex analog +// by replacing all entries with a 2-vector that has zero imaginary part. +function complex(list) = + is_num(list) ? [list,0] : + [for(entry=list) is_num(entry) ? [entry,0] : complex(entry)]; + + // Function: c_mul() // Usage: // c = c_mul(z1,z2) // Description: -// Multiplies two complex numbers represented by 2D vectors. -// Returns a complex number as a 2D vector [REAL, IMAGINARY]. +// Multiplies two complex numbers, vectors or matrices, where complex numbers +// or entries are represented as vectors: [REAL, IMAGINARY]. Note that all +// entries in both arguments must be complex. // Arguments: -// z1 = First complex number, given as a 2D vector [REAL, IMAGINARY] -// z2 = Second complex number, given as a 2D vector [REAL, IMAGINARY] -function c_mul(z1,z2) = - assert( is_matrix([z1,z2],2,2), "Complex numbers should be represented by 2D vectors" ) +// z1 = First complex number, vector or matrix +// z2 = Second complex number, vector or matrix + +function _split_complex(data) = + is_vector(data,2) ? data + : is_num(data[0][0]) ? [data*[1,0], data*[0,1]] + : [ + [for(vec=data) vec * [1,0]], + [for(vec=data) vec * [0,1]] + ]; + +function _combine_complex(data) = + is_vector(data,2) ? data + : is_num(data[0][0]) ? [for(i=[0:len(data[0])-1]) [data[0][i],data[1][i]]] + : [for(i=[0:1:len(data[0])-1]) + [for(j=[0:1:len(data[0][0])-1]) + [data[0][i][j], data[1][i][j]]]]; + +function _c_mul(z1,z2) = [ z1.x*z2.x - z1.y*z2.y, z1.x*z2.y + z1.y*z2.x ]; +function c_mul(z1,z2) = + is_matrix([z1,z2],2,2) ? _c_mul(z1,z2) : + _combine_complex(_c_mul(_split_complex(z1), _split_complex(z2))); + + // Function: c_div() // Usage: // x = c_div(z1,z2) @@ -1493,7 +1526,52 @@ function c_div(z1,z2) = let(den = z2.x*z2.x + z2.y*z2.y) [(z1.x*z2.x + z1.y*z2.y)/den, (z1.y*z2.x - z1.x*z2.y)/den]; -// For the sake of consistence with Q_mul and vmul, c_mul should be called C_mul + +// Function: c_conj() +// Usage: +// w = c_conj(z) +// Description: +// Computes the complex conjugate of the input, which can be a complex number, +// complex vector or complex matrix. +function c_conj(z) = + is_vector(z,2) ? [z.x,-z.y] : + [for(entry=z) c_conj(entry)]; + +// Function: c_real() +// Usage: +// x = c_real(z) +// Description: +// Returns real part of a complex number, vector or matrix. +function c_real(z) = + is_vector(z,2) ? z.x + : is_num(z[0][0]) ? z*[1,0] + : [for(vec=z) vec * [1,0]]; + +// Function: c_imag() +// Usage: +// x = c_imag(z) +// Description: +// Returns imaginary part of a complex number, vector or matrix. +function c_imag(z) = + is_vector(z,2) ? z.y + : is_num(z[0][0]) ? z*[0,1] + : [for(vec=z) vec * [0,1]]; + + +// Function: c_ident() +// Usage: +// I = c_ident(n) +// Description: +// Produce an n by n complex identity matrix +function c_ident(n) = [for (i = [0:1:n-1]) [for (j = [0:1:n-1]) (i==j)?[1,0]:[0,0]]]; + +// Function: c_norm() +// Usage: +// n = c_norm(z) +// Description: +// Compute the norm of a complex number or vector. +function c_norm(z) = norm_fro(z); + // Section: Polynomials @@ -1539,12 +1617,12 @@ function quadratic_roots(a,b,c,real=false) = // where a_n is the z^n coefficient. Polynomial coefficients are real. // The result is a number if `z` is a number and a complex number otherwise. function polynomial(p,z,k,total) = -    is_undef(k) -    ?   assert( is_vector(p) , "Input polynomial coefficients must be a vector." ) -        assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." ) -        polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0]) -    : k==len(p) ? total -    : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : c_mul(total,z)+[p[k],0]); + is_undef(k) + ? assert( is_vector(p) , "Input polynomial coefficients must be a vector." ) + assert( is_finite(z) || is_vector(z,2), "The value of `z` must be a real or a complex number." ) + polynomial( _poly_trim(p), z, 0, is_num(z) ? 0 : [0,0]) + : k==len(p) ? total + : polynomial(p,z,k+1, is_num(z) ? total*z+p[k] : c_mul(total,z)+[p[k],0]); // Function: poly_mult() // Usage: @@ -1554,12 +1632,12 @@ function polynomial(p,z,k,total) = // Given a list of polynomials represented as real coefficient lists, with the highest degree coefficient first, // computes the coefficient list of the product polynomial. function poly_mult(p,q) = -    is_undef(q) ? -        len(p)==2 + is_undef(q) ? + len(p)==2 ? poly_mult(p[0],p[1]) -        : poly_mult(p[0], poly_mult(select(p,1,-1))) -    : -    assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult") + : poly_mult(p[0], poly_mult(select(p,1,-1))) + : + assert( is_vector(p) && is_vector(q),"Invalid arguments to poly_mult") p*p==0 || q*q==0 ? [0] : _poly_trim(convolve(p,q)); diff --git a/skin.scad b/skin.scad index de15cda..92536aa 100644 --- a/skin.scad +++ b/skin.scad @@ -407,7 +407,7 @@ function skin(profiles, slices, refine=1, method="direct", sampling, caps, close assert(methodlistok==[], str("method list contains invalid method at ",methodlistok)) assert(len(method) == profcount,"Method list is the wrong length") assert(in_list(sampling,["length","segment"]), "sampling must be set to \"length\" or \"segment\"") - assert(sampling=="segment" || (!in_list("distance",method) && !in_list("fast_distance") && !in_list("tangent",method)), "sampling is set to \"length\" which is only allowed with methods \"direct\" and \"reindex\"") + assert(sampling=="segment" || (!in_list("distance",method) && !in_list("fast_distance",method) && !in_list("tangent",method)), "sampling is set to \"length\" which is only allowed with methods \"direct\" and \"reindex\"") assert(capsOK, "caps must be boolean or a list of two booleans") assert(!closed || !caps, "Cannot make closed shape with caps") let( diff --git a/tests/test_math.scad b/tests/test_math.scad index b0e1782..ae9801d 100644 --- a/tests/test_math.scad +++ b/tests/test_math.scad @@ -824,19 +824,78 @@ module test_lcm() { test_lcm(); -module test_C_times() { - assert_equal(C_times([4,5],[9,-4]), [56,29]); - assert_equal(C_times([-7,2],[24,3]), [-174, 27]); +module test_c_mul() { + assert_equal(c_mul([4,5],[9,-4]), [56,29]); + assert_equal(c_mul([-7,2],[24,3]), [-174, 27]); + assert_equal(c_mul([3,4], [[3,-7], [4,9], [4,8]]), [[37,-9],[-24,43], [-20,40]]); + assert_equal(c_mul([[3,-7], [4,9], [4,8]], [[1,1],[3,4],[-3,4]]), [-58,31]); + M = [ + [ [3,4], [9,-1], [4,3] ], + [ [2,9], [4,9], [3,-1] ] + ]; + assert_equal(c_mul(M, [ [3,4], [4,4],[5,5]]), [[38,91], [-30, 97]]); + assert_equal(c_mul([[4,4],[9,1]], M), [[5,111],[67,117], [32,22]]); + assert_equal(c_mul(M,transpose(M)), [ [[80,30], [30, 117]], [[30,117], [-134, 102]]]); + assert_equal(c_mul(transpose(M),M), [ [[-84,60],[-42,87],[15,50]], [[-42,87],[15,54],[60,46]], [[15,50],[60,46],[15,18]]]); } -test_C_times(); +test_c_mul(); -module test_C_div() { - assert_equal(C_div([56,29],[9,-4]), [4,5]); - assert_equal(C_div([-174,27],[-7,2]), [24,3]); +module test_c_div() { + assert_equal(c_div([56,29],[9,-4]), [4,5]); + assert_equal(c_div([-174,27],[-7,2]), [24,3]); } -test_C_div(); +test_c_div(); +module test_c_conj(){ + assert_equal(c_conj([3,4]), [3,-4]); + assert_equal(c_conj( [ [2,9], [4,9], [3,-1] ]), [ [2,-9], [4,-9], [3,1] ]); + M = [ + [ [3,4], [9,-1], [4,3] ], + [ [2,9], [4,9], [3,-1] ] + ]; + Mc = [ + [ [3,-4], [9,1], [4,-3] ], + [ [2,-9], [4,-9], [3,1] ] + ]; + assert_equal(c_conj(M), Mc); +} +test_c_conj(); + +module test_c_real(){ + M = [ + [ [3,4], [9,-1], [4,3] ], + [ [2,9], [4,9], [3,-1] ] + ]; + assert_equal(c_real(M), [[3,9,4],[2,4,3]]); + assert_equal(c_real( [ [3,4], [9,-1], [4,3] ]), [3,9,4]); + assert_equal(c_real([3,4]),3); +} +test_c_real(); + + +module test_c_imag(){ + M = [ + [ [3,4], [9,-1], [4,3] ], + [ [2,9], [4,9], [3,-1] ] + ]; + assert_equal(c_imag(M), [[4,-1,3],[9,9,-1]]); + assert_equal(c_imag( [ [3,4], [9,-1], [4,3] ]), [4,-1,3]); + assert_equal(c_imag([3,4]),4); +} +test_c_imag(); + + +module test_c_ident(){ + assert_equal(c_ident(3), [[[1, 0], [0, 0], [0, 0]], [[0, 0], [1, 0], [0, 0]], [[0, 0], [0, 0], [1, 0]]]); +} +test_c_ident(); + +module test_c_norm(){ + assert_equal(c_norm([3,4]), 5); + assert_approx(c_norm([[3,4],[5,6]]), 9.273618495495704); +} +test_c_norm(); module test_back_substitute(){ R = [[12,4,3,2],