Whilst working through issues with porting Sage to FreeBSD, I've run into a lack of complex trig/exponential functions - which Sage (specifically sage-x.x.x/sage/ext/interpreters/wrapper_cdf.pyx) expects.
BTW, how is wrapper_cdf.pyx generated? There's an 'auto-generated' header but I can't find anything obvious that generates it. How much does sage care about the correct handling of exceptional IEEE conditions in libm functions? Things like correctly propogating NaN, infinity and the sign of zero. The easy way to implement the missing-in-FreeBSD functions is to directly implement mathematical equivalences like: cexp(z) = z.r * (cos(z.i) + I * sin(z.i)) csin(z) = sin(z.r) * cos(I * z.i) + cos(z.r) * sin(I * z.i) csinh(z) = sinh(z.r) * cos(z.i) + I * cosh(z.r) * sin(z.i) ctan(z) = (tan(z.r) + tan(I * z.i)) / (1 - tan(z.r) * tan(I * z.i)) clog(z) = log(sqrt(z.r*z.r + z.i*z.i)) + I * atan(z.i/z.r) cpow(x, y) = cexp(y * clog(x)) but these do not necessarily behave well for exceptional numbers and the latter 3 can lose significant amounts of precision and/or suffer from unexpected underflow/overflow. The options I can think of are: 1) Chop the unimplemented functions out of wrapper_cdf (and related bits of sage). 2) Create a simplistic implementation (as per the above) 3) Create or borrow a "proper" implementation that does the appropriate exception checking and semi-numerical decomposition. The latter two options additionally offer the choice of directly embedding the code into wrapper_cdf.so (probably the better approach for option 2) or creating a new FreeBSD-only spkg that creates a .so that wrapper_cdf.so needs to depend on (probably the better approach for option 3). Options for "borrowing" code include NetBSD and glibc - though neither do a complete job (eg the NetBSD code appears to skimp on the checks on input values and glibc-2.9 uses "cpow(x, y) = cexp(y * clog(x))"). Note that doing it "correctly" is non-trivial - relevant decomposition, checks and handling zero, NaN and Inf typically add a page of C code. Even something as simple as cabs(z) aka hypot(z) needs about 70 lines of C to do properly. -- Peter Jeremy
pgp5Hu6J03EJz.pgp
Description: PGP signature
