|
Hi Fred, could you please run the attached debug version of ComplexCell.cc? I am getting: 5J3 | 4J¯1 ⎕CT is: 1e-13 modulus (A) is: (5,3) A=0 is: (0,0) A+A=0 is: (5,3) B÷A+A=0 is: (0.5,-0.5) ⌊B÷A+A=0 is: (1,-1) A×⌊B÷A+A=0 is: (8,-2) B-A×⌊B÷A+A=0 is: (-4,1) ¯4J1 The above printout are the intermediate results of the following computation in bif_resudue(): // // op Z before Z after new (Z) IntCell(0); // Z←0 any 0 Z->bif_equal(Z, A); // Z←A=Z 0 A=0 Z->bif_add(Z, A); // Z←A+Z A=0 A+A=0 Z->bif_divide(Z, this); // Z←B÷Z A+A=0 B÷A+A=0 Z->bif_floor(Z); // Z←⌊Z B÷A+A=0 ⌊B÷A+A=0 Z->bif_multiply(Z, A); // Z←A×Z ⌊B÷A+A=0 A×⌊B÷A+A=0 Z->bif_subtract(Z, this); // Z←A-Z A×⌊B÷A+A=0 B-A×⌊B÷A+A=0 which I believe is a faithful (apart from errors) implementation of Z←B-A×⌊B÷A+A=0 as given in the McDonnell paper. The key issue are probably the following out put lines: B÷A+A=0 is: (0.5,-0.5) ⌊B÷A+A=0 is: (1,-1) Best Regards, /// Jürgen On 06/24/2017 01:48 AM, Frederick Pitts
wrote:
|
/*
This file is part of GNU APL, a free implementation of the
ISO/IEC Standard 13751, "Programming Language APL, Extended"
Copyright (C) 2008-2016 Dr. Jürgen Sauermann
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <errno.h>
#include <fenv.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "ComplexCell.hh"
#include "ErrorCode.hh"
#include "FloatCell.hh"
#include "IntCell.hh"
#include "Output.hh"
#include "Workspace.hh"
#include "Cell.icc"
//-----------------------------------------------------------------------------
ComplexCell::ComplexCell(APL_Complex c)
{
value.cval_r = c.real();
value2.cval_i = c.imag();
}
//-----------------------------------------------------------------------------
ComplexCell::ComplexCell(APL_Float r, APL_Float i)
{
value.cval_r = r;
value2.cval_i = i;
}
//-----------------------------------------------------------------------------
APL_Float
ComplexCell::get_real_value() const
{
return value.cval_r;
}
//-----------------------------------------------------------------------------
APL_Float
ComplexCell::get_imag_value() const
{
return value2.cval_i;
}
//-----------------------------------------------------------------------------
bool
ComplexCell::is_near_int() const
{
return Cell::is_near_int(value.cval_r) &&
Cell::is_near_int(value2.cval_i);
}
//-----------------------------------------------------------------------------
bool
ComplexCell::is_near_zero() const
{
if (value.cval_r >= INTEGER_TOLERANCE) return false;
if (value.cval_r <= -INTEGER_TOLERANCE) return false;
if (value2.cval_i >= INTEGER_TOLERANCE) return false;
if (value2.cval_i <= -INTEGER_TOLERANCE) return false;
return true;
}
//-----------------------------------------------------------------------------
bool
ComplexCell::is_near_one() const
{
if (value.cval_r > (1.0 + INTEGER_TOLERANCE)) return false;
if (value.cval_r < (1.0 - INTEGER_TOLERANCE)) return false;
if (value2.cval_i > INTEGER_TOLERANCE) return false;
if (value2.cval_i < -INTEGER_TOLERANCE) return false;
return true;
}
//-----------------------------------------------------------------------------
bool
ComplexCell::is_near_real() const
{
const APL_Float B = REAL_TOLERANCE*REAL_TOLERANCE;
const double I = value2.cval_i * value2.cval_i;
if (I < B) return true; // I is absolutely small
const double R = value.cval_r * value.cval_r;
if (I < R*B) return true; // I is relatively small
return false;
}
//-----------------------------------------------------------------------------
bool
ComplexCell::equal(const Cell & A, APL_Float qct) const
{
if (!A.is_numeric()) return false;
return tolerantly_equal(A.get_complex_value(), get_complex_value(), qct);
}
//-----------------------------------------------------------------------------
// monadic build-in functions...
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_factorial(Cell * Z) const
{
if (is_near_real())
{
const FloatCell fc(get_real_value());
fc.bif_factorial(Z);
return E_NO_ERROR;
}
const APL_Float zr = get_real_value() + 1.0;
const APL_Float zi = get_imag_value();
const APL_Complex z(zr, zi);
new (Z) ComplexCell(gamma(get_real_value(), get_imag_value()));
if (errno) return E_DOMAIN_ERROR;
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_conjugate(Cell * Z) const
{
new (Z) ComplexCell(conj(cval()));
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_negative(Cell * Z) const
{
new (Z) ComplexCell(-cval());
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_direction(Cell * Z) const
{
const APL_Float mag = abs(cval());
if (mag == 0.0) return IntCell::zv(Z, 0);
else return ComplexCell::zv(Z, cval()/mag);
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_magnitude(Cell * Z) const
{
return FloatCell::zv(Z, abs(cval()));
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_reciprocal(Cell * Z) const
{
if (is_near_zero()) return E_DOMAIN_ERROR;
if (is_near_real())
{
const APL_Float z = 1.0/value.cval_r;
if (!isfinite(z)) return E_DOMAIN_ERROR;
return FloatCell::zv(Z, z);
}
const APL_Complex z = 1.0 / cval();
if (!isfinite(z.real())) return E_DOMAIN_ERROR;
if (!isfinite(z.imag())) return E_DOMAIN_ERROR;
return zv(Z, z);
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_roll(Cell * Z) const
{
const APL_Integer qio = Workspace::get_IO();
if (!is_near_int()) return E_DOMAIN_ERROR;
const APL_Integer set_size = get_checked_near_int();
if (set_size <= 0) return E_DOMAIN_ERROR;
const uint64_t rnd = Workspace::get_RL(set_size);
return IntCell::zv(Z, qio + (rnd % set_size));
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_pi_times(Cell * Z) const
{
new (Z) ComplexCell(M_PI * cval());
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_pi_times_inverse(Cell * Z) const
{
new (Z) ComplexCell(cval() / M_PI);
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_ceiling(Cell * Z) const
{
const APL_Float cr = ceil(value.cval_r); // cr ⥠value.cval_r
const APL_Float Dr = cr - value.cval_r; // Dr ⥠0
const APL_Float ci = ceil(value2.cval_i); // ci ⥠value2.cval_i
const APL_Float Di = ci - value2.cval_i; // Di ⥠0
const APL_Float D = Dr + Di; // 0 ⤠D < 2
// ISO: if D is tolerantly less than 1 return fr + 0J1Ãfi
// IBM: if D is less than 1 return fr + 0J1Ãfi
// However, ISO examples follow IBM (and so do we)
//
// if (D < 1.0 + Workspace::get_CT()) return zv(Z, cr, ci); // ISO
if (D < 1.0) return zv(Z, cr, ci); // IBM and examples in ISO
if (Di > Dr) return zv(Z, cr, ci - 1.0);
else return zv(Z, cr - 1.0, ci);
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_floor(Cell * Z) const
{
const APL_Float fr = floor(value.cval_r); // fr ⤠value.cval_r
const APL_Float Dr = value.cval_r - fr; // 0 ⤠Dr < 1
const APL_Float fi = floor(value2.cval_i); // fi ⤠value2.cval_i
const APL_Float Di = value2.cval_i - fi; // 0 ⤠Di < 1
const APL_Float D = Dr + Di; // 0 ⤠D < 2
if (D < 1.0 - Workspace::get_CT()) return zv(Z, fr, fi); // ISO
if (Dr < Di) return zv(Z, fr, fi + 1.0);
else return zv(Z, fr + 1.0, fi);
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_exponential(Cell * Z) const
{
new (Z) ComplexCell(exp(cval()));
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_nat_log(Cell * Z) const
{
new (Z) ComplexCell(log(cval()));
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
// dyadic build-in functions...
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_add(Cell * Z, const Cell * A) const
{
new (Z) ComplexCell(A->get_complex_value() + get_complex_value());
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_subtract(Cell * Z, const Cell * A) const
{
new (Z) ComplexCell(A->get_complex_value() - get_complex_value());
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_multiply(Cell * Z, const Cell * A) const
{
const APL_Complex z = A->get_complex_value() * cval();
if (!isfinite(z.real())) return E_DOMAIN_ERROR;
if (!isfinite(z.imag())) return E_DOMAIN_ERROR;
new (Z) ComplexCell(z);
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_divide(Cell * Z, const Cell * A) const
{
if (cval().real() == 0.0 && cval().imag() == 0.0) // A ÷ 0
{
if (A->get_real_value() != 0.0) return E_DOMAIN_ERROR;
if (A->get_imag_value() != 0.0) return E_DOMAIN_ERROR;
return IntCell::zv(Z, 1); // 0÷0 is 1 in APL
}
new (Z) ComplexCell(A->get_complex_value() / get_complex_value());
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_residue(Cell * Z, const Cell * A) const
{
if (!A->is_numeric()) return E_DOMAIN_ERROR;
const APL_Complex mod = A->get_complex_value();
const APL_Complex b = get_complex_value();
// if A is zero , return B
//
if (mod.real() == 0.0 && mod.imag() == 0.0)
return zv(Z, b);
// IBM: if B is zero , return 0
//
if (b.real() == 0.0 && b.imag() == 0.0)
return IntCell::z0(Z);
// compliant implementation: B-AÃâB÷A+A=0
// The b=0 case may still exist due to âCT
//
// // op Z before Z after
CERR << std::setprecision(14);
CERR << "âCT is: " << Workspace::get_CT() << endl;
CERR << "modulus (A) is: " << A->get_complex_value() << endl;
new (Z) IntCell(0); // Zâ0 any 0
Z->bif_equal(Z, A); // ZâA=Z 0 A=0
CERR << "A=0 is: " << Z->get_complex_value() << endl;
Z->bif_add(Z, A); // ZâA+Z A=0 A+A=0
CERR << "A+A=0 is: " << Z->get_complex_value() << endl;
Z->bif_divide(Z, this); // ZâB÷Z A+A=0 B÷A+A=0
CERR << "B÷A+A=0 is: " << Z->get_complex_value() << endl;
Z->bif_floor(Z); // ZââZ B÷A+A=0 âB÷A+A=0
CERR << "âB÷A+A=0 is: " << Z->get_complex_value() << endl;
Z->bif_multiply(Z, A); // ZâAÃZ âB÷A+A=0 AÃâB÷A+A=0
CERR << "AÃâB÷A+A=0 is: " << Z->get_complex_value() << endl;
Z->bif_subtract(Z, this); // ZâA-Z AÃâB÷A+A=0 B-AÃâB÷A+A=0
CERR << "B-AÃâB÷A+A=0 is: " << Z->get_complex_value() << endl;
return E_NO_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_maximum(Cell * Z, const Cell * A) const
{
// maximum of complex numbers gives DOMAN ERROR if one of the cells
// is not near-real.
//
if (!is_near_real()) return E_DOMAIN_ERROR;
if (!A->is_near_real()) return E_DOMAIN_ERROR;
const APL_Float breal = get_real_value();
if (A->is_integer_cell())
{
const APL_Integer a = A->get_int_value();
return a >= breal ? IntCell::zv(Z, a) : FloatCell::zv(Z, breal);
}
// both A and B are near-real, therefore they must be numeric and have no
// non-0 imag part
//
const APL_Float a = A->get_real_value();
return FloatCell::zv(Z, a >= breal ? a : breal);
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_minimum(Cell * Z, const Cell * A) const
{
// minimum of complex numbers gives DOMAN ERROR if one of the cells
// is not near real.
//
if (!is_near_real()) return E_DOMAIN_ERROR;
if (!A->is_near_real()) return E_DOMAIN_ERROR;
const APL_Float breal = get_real_value();
if (A->is_integer_cell())
{
const APL_Integer a = A->get_int_value();
return a <= breal ? IntCell::zv(Z, a) : FloatCell::zv(Z, breal);
}
// both A and B are near-real, therefore they must be numeric and have no
// non-0 imag part
//
const APL_Float a = A->get_real_value();
return FloatCell::zv(Z, a <= breal ? a : breal);
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_equal(Cell * Z, const Cell * A) const
{
const APL_Float qct = Workspace::get_CT();
if (A->is_complex_cell())
{
const APL_Complex diff = A->get_complex_value() - get_complex_value();
if (diff.real() > qct) return IntCell::zv(Z, 0);
if (diff.real() < -qct) return IntCell::zv(Z, 0);
if (diff.imag() > qct) return IntCell::zv(Z, 0);
if (diff.imag() < -qct) return IntCell::zv(Z, 0);
return IntCell::zv(Z, 1);
}
if (A->is_numeric())
{
if (get_imag_value() > qct) return IntCell::zv(Z, 0);
if (get_imag_value() < -qct) return IntCell::zv(Z, 0);
const APL_Float diff = A->get_real_value() - get_real_value();
if (diff > qct) return IntCell::zv(Z, 0);
if (diff < -qct) return IntCell::zv(Z, 0);
return IntCell::zv(Z, 1);
}
return IntCell::zv(Z, 1);
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_power(Cell * Z, const Cell * A) const
{
// some A to the complex B-th power
//
if (!A->is_numeric()) return E_DOMAIN_ERROR;
const APL_Float ar = A->get_real_value();
const APL_Float ai = A->get_imag_value();
// 1. A == 0
//
if (ar == 0.0 && ai == 0.0)
{
if (cval().real() == 0.0) return IntCell::z1(Z); // 0â0 is 1
if (cval().imag() > 0) return IntCell::z0(Z); // 0âN is 0
return E_DOMAIN_ERROR; // 0â¯N = 1÷0
}
// 2. complex result
{
const APL_Complex a(ar, ai);
const APL_Complex z = pow(a, cval());
if (!isfinite(z.real())) return E_DOMAIN_ERROR;
if (!isfinite(z.imag())) return E_DOMAIN_ERROR;
new (Z) ComplexCell(z);
return E_NO_ERROR;
}
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_logarithm(Cell * Z, const Cell * A) const
{
// ISO p. 88
//
if (!A->is_numeric()) return E_DOMAIN_ERROR;
if (get_real_value() == A->get_real_value() &&
get_imag_value() == A->get_imag_value()) return IntCell::z1(Z);
if (value.cval_r == 0.0 && value2.cval_i == 0.0) return E_DOMAIN_ERROR;
if (A->is_near_one()) return E_DOMAIN_ERROR;
if (A->is_real_cell())
{
const APL_Complex z = log(cval()) / log(A->get_real_value());
if (!isfinite(z.real())) return E_DOMAIN_ERROR;
if (!isfinite(z.imag())) return E_DOMAIN_ERROR;
return zv(Z, z);
}
if (A->is_complex_cell())
{
const APL_Complex z = log(cval()) / log(A->get_complex_value());
if (!isfinite(z.real())) return E_DOMAIN_ERROR;
if (!isfinite(z.imag())) return E_DOMAIN_ERROR;
return zv(Z, z);
}
return E_DOMAIN_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_circle_fun(Cell * Z, const Cell * A) const
{
if (!A->is_near_int()) return E_DOMAIN_ERROR;
const APL_Integer fun = A->get_checked_near_int();
new (Z) FloatCell(0); // prepare for DOMAIN ERROR
const ErrorCode ret = do_bif_circle_fun(Z, fun, cval());
if (!Z->is_finite()) return E_DOMAIN_ERROR;
return ret;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::bif_circle_fun_inverse(Cell * Z, const Cell * A) const
{
if (!A->is_near_int()) return E_DOMAIN_ERROR;
const APL_Integer fun = A->get_checked_near_int();
new (Z) FloatCell(0); // prepare for DOMAIN ERROR
ErrorCode ret = E_DOMAIN_ERROR;
switch(fun)
{
case 1: case -1:
case 2: case -2:
case 3: case -3:
case 4: case -4:
case 5: case -5:
case 6: case -6:
case 7: case -7:
ret = do_bif_circle_fun(Z, -fun, cval());
if (!Z->is_finite()) return E_DOMAIN_ERROR;
return ret;
case -10: // +A (conjugate) is self-inverse
ret = do_bif_circle_fun(Z, fun, cval());
if (!Z->is_finite()) return E_DOMAIN_ERROR;
return ret;
default: return E_DOMAIN_ERROR;
}
// not reached
return E_DOMAIN_ERROR;
}
//-----------------------------------------------------------------------------
ErrorCode
ComplexCell::do_bif_circle_fun(Cell * Z, int fun, APL_Complex b)
{
switch(fun)
{
case -12: {
ComplexCell cb(-b.imag(), b.real());
cb.bif_exponential(Z);
}
return E_NO_ERROR;
case -11: return zv(Z, -b.imag(), b.real());
case -10: return zv(Z, b.real(), -b.imag());
case -9: return zv(Z, b );
case -7: // arctanh(z) = 0.5 (ln(1 + z) - ln(1 - z))
{
const APL_Complex b1 = ONE() + b;
const APL_Complex b_1 = ONE() - b;
const APL_Complex log_b1 = log(b1);
const APL_Complex log_b_1 = log(b_1);
const APL_Complex diff = log_b1 - log_b_1;
const APL_Complex half = 0.5 * diff;
return zv(Z, half);
}
case -6: // arccosh(z) = ln(z + sqrt(z + 1) sqrt(z - 1))
{
const APL_Complex b1 = b + ONE();
const APL_Complex b_1 = b - ONE();
const APL_Complex root1 = sqrt(b1);
const APL_Complex root_1 = sqrt(b_1);
const APL_Complex prod = root1 * root_1;
const APL_Complex sum = b + prod;
const APL_Complex loga = log(sum);
return zv(Z, loga);
}
case -5: // arcsinh(z) = ln(z + sqrt(z^2 + 1))
{
const APL_Complex b2 = b*b;
const APL_Complex b2_1 = b2 + ONE();
const APL_Complex root = sqrt(b2_1);
const APL_Complex sum = b + root;
const APL_Complex loga = log(sum);
return zv(Z, loga);
}
case -4: if (b.real() >= 0 ||
(b.real() > -1 && Cell::is_near_zero(b.imag()))
) return zv(Z, sqrt(b*b - 1.0));
else return zv(Z, -sqrt(b*b - 1.0));
case -3: // arctan(z) = i/2 (ln(1 - iz) - ln(1 + iz))
{
const APL_Complex iz = APL_Complex(- b.imag(), b.real());
const APL_Complex piz = ONE() + iz;
const APL_Complex niz = ONE() - iz;
const APL_Complex log_piz = log(piz);
const APL_Complex log_niz = log(niz);
const APL_Complex diff = log_niz - log_piz;
const APL_Complex prod = APL_Complex(0, 0.5) * diff;
return zv(Z, prod);
}
case -2: // arccos(z) = -i (ln( z + sqrt(z^2 - 1)))
{
const APL_Complex b2 = b*b;
const APL_Complex diff = b2 - ONE();
const APL_Complex root = sqrt(diff);
const APL_Complex sum = b + root;
const APL_Complex loga = log(sum);
const APL_Complex prod = MINUS_i() * loga;
return zv(Z, prod);
}
case -1: // arcsin(z) = -i (ln(iz + sqrt(1 - z^2)))
{
const APL_Complex b2 = b*b;
const APL_Complex diff = ONE() - b2;
const APL_Complex root = sqrt(diff);
const APL_Complex sum = APL_Complex(-b.imag(), b.real())
+ root;
const APL_Complex loga = log(sum);
const APL_Complex prod = MINUS_i() * loga;
return zv(Z, prod);
}
case 0: return zv(Z, sqrt(1.0 - b*b));
case 1: return zv(Z, sin (b));
case 2: return zv(Z, cos (b));
case 3: return zv(Z, tan (b));
case 4: return zv(Z, sqrt(1.0 + b*b));
case 5: return zv(Z, sinh (b));
case 6: return zv(Z, cosh (b));
case 7: return zv(Z, tanh (b));
case -8: // ¯8âB is - 8âB
case 8: {
// Let b = X + iY and sq = (¯1 - R*2)*.5
// Then 8âX+JY: Z is: + sq if X > 0 and Y > 0
// or X = 0 and Y > 1
// or X < 0 and Y ⥠0, and
// - sq otherwise
bool pos_8 = false;
if (b.real() > 0) { if (b.imag() > 0) pos_8 = true; }
else if (b.real() == 0) { if (b.imag() > 1) pos_8 = true; }
else { if (b.imag() >= 0) pos_8 = true; }
if (fun == 8) pos_8 = ! pos_8;
const APL_Complex sq = sqrt(-(1.0 + b*b)); // (¯1 - R*2)*.5
return zv(Z, pos_8 ? sq : -sq);
}
case 9: return FloatCell::zv(Z, b.real());
case 10: return FloatCell::zv(Z, sqrt(b.real()*b.real()
+ b.imag()* b.imag()));
case 11: return FloatCell::zv(Z, b.imag());
case 12: return FloatCell::zv(Z, atan2(b.imag(), b.real()));
}
// invalid fun
//
return E_DOMAIN_ERROR;
}
//-----------------------------------------------------------------------------
APL_Complex
ComplexCell::gamma(APL_Float x, const APL_Float y)
{
if (x < 0.5)
return M_PI / sin(M_PI * APL_Complex(x, y)) * gamma(1 - x, -y);
// coefficients for lanczos approximation of the gamma function.
//
#define p0 APL_Complex( 1.000000000190015 )
#define p1 APL_Complex( 76.18009172947146 )
#define p2 APL_Complex(-86.50532032941677 )
#define p3 APL_Complex( 24.01409824083091 )
#define p4 APL_Complex( -1.231739572450155 )
#define p5 APL_Complex( 1.208650973866179E-3)
#define p6 APL_Complex( -5.395239384953E-6 )
errno = 0;
feclearexcept(FE_ALL_EXCEPT);
x += 1.0;
const APL_Complex z(x, y);
const APL_Complex ret( (sqrt(2*M_PI) / z)
* (p0 +
p1 / APL_Complex(x + 1.0, y) +
p2 / APL_Complex(x + 2.0, y) +
p3 / APL_Complex(x + 3.0, y) +
p4 / APL_Complex(x + 4.0, y) +
p5 / APL_Complex(x + 5.0, y) +
p6 / APL_Complex(x + 6.0, y))
* pow(z + 5.5, z + 0.5)
* exp(-(z + 5.5))
);
// errno may be set and is checked by the caller
return ret;
#undef p0
#undef p1
#undef p2
#undef p3
#undef p4
#undef p5
#undef p6
}
//=============================================================================
// throw/nothrow boundary. Functions above MUST NOT (directly or indirectly)
// throw while funcions below MAY throw.
//=============================================================================
#include "Error.hh" // throws
//-----------------------------------------------------------------------------
bool
ComplexCell::get_near_bool() const
{
if (value2.cval_i > INTEGER_TOLERANCE) DOMAIN_ERROR;
if (value2.cval_i < -INTEGER_TOLERANCE) DOMAIN_ERROR;
if (value.cval_r > INTEGER_TOLERANCE) // 1 or invalid
{
if (value.cval_r < (1.0 - INTEGER_TOLERANCE)) DOMAIN_ERROR;
if (value.cval_r > (1.0 + INTEGER_TOLERANCE)) DOMAIN_ERROR;
return true;
}
else
{
if (value.cval_r < -INTEGER_TOLERANCE) DOMAIN_ERROR;
return false;
}
}
//-----------------------------------------------------------------------------
APL_Integer
ComplexCell::get_near_int() const
{
// if (value2.cval_i > qct) DOMAIN_ERROR;
// if (value2.cval_i < -qct) DOMAIN_ERROR;
const double val = value.cval_r;
const double result = round(val);
const double diff = val - result;
if (diff > INTEGER_TOLERANCE) DOMAIN_ERROR;
if (diff < -INTEGER_TOLERANCE) DOMAIN_ERROR;
return APL_Integer(result + 0.3);
}
//-----------------------------------------------------------------------------
bool
ComplexCell::greater(const Cell & other) const
{
switch(other.get_cell_type())
{
case CT_CHAR: return true;
case CT_INT:
case CT_FLOAT:
case CT_COMPLEX: break;
case CT_POINTER: return false;
case CT_CELLREF: DOMAIN_ERROR;
default: Assert(0 && "Bad celltype");
}
const Comp_result comp = compare(other);
if (comp == COMP_EQ) return this > &other;
return (comp == COMP_GT);
}
//-----------------------------------------------------------------------------
Comp_result
ComplexCell::compare(const Cell & other) const
{
if (other.is_character_cell()) return COMP_GT;
if (!other.is_numeric()) DOMAIN_ERROR;
const APL_Float qct = Workspace::get_CT();
if (equal(other, qct)) return COMP_EQ;
APL_Float areal = other.get_real_value();
APL_Float breal = get_real_value();
// we compare complex numbers by their real part unless the real parsts
// are equal. In that case we compare the imag parts. The reason for
// comparing this way is compatibility with the comparison of real numbers
//
if (tolerantly_equal(areal, breal, qct))
{
areal = other.get_imag_value();
breal = get_imag_value();
}
return (breal < areal) ? COMP_LT : COMP_GT;
}
//-----------------------------------------------------------------------------
PrintBuffer
ComplexCell::character_representation(const PrintContext & pctx) const
{
if (pctx.get_PP() < MAX_Quad_PP)
{
// 10âget_PP()
//
APL_Float ten_to_PP = 1.0;
loop(p, pctx.get_PP()) ten_to_PP *= 10.0;
// lrm p. 13: In J notation, the real or imaginary part is not
// displayed if it is less than the other by more than âPP orders
// of magnitude (unless âPP is at its maximum).
//
const APL_Float pos_real = value.cval_r < 0
? -value.cval_r : value.cval_r;
const APL_Float pos_imag = value2.cval_i < 0
? -value2.cval_i : value2.cval_i;
if (pos_real >= pos_imag) // pos_real dominates pos_imag
{
if (pos_real > pos_imag*ten_to_PP)
{
const FloatCell real_cell(value.cval_r);
return real_cell.character_representation(pctx);
}
}
else // pos_imag dominates pos_real
{
if (pos_imag > pos_real*ten_to_PP)
{
const FloatCell imag_cell(value2.cval_i);
PrintBuffer ret = imag_cell.character_representation(pctx);
ret.pad_l(UNI_ASCII_J, 1);
ret.pad_l(UNI_ASCII_0, 1);
ret.get_info().flags |= CT_COMPLEX;
ret.get_info().imag_len = 1 + ret.get_info().real_len;
ret.get_info().int_len = 1;
ret.get_info().fract_len = 0;
ret.get_info().real_len = 1;
return ret;
}
}
}
bool scaled_real = pctx.get_scaled(); // may be changed by print function
UCS_string ucs(value.cval_r, scaled_real, pctx);
ColInfo info;
info.flags |= CT_COMPLEX;
if (scaled_real) info.flags |= real_has_E;
int int_fract = ucs.size();
info.real_len = ucs.size();
info.int_len = ucs.size();
loop(u, ucs.size())
{
if (ucs[u] == UNI_ASCII_FULLSTOP)
{
info.int_len = u;
if (!scaled_real) break;
continue;
}
if (ucs[u] == UNI_ASCII_E)
{
if (info.int_len > u) info.int_len = u;
int_fract = u;
break;
}
}
info.fract_len = int_fract - info.int_len;
if (!is_near_real())
{
ucs.append(UNI_ASCII_J);
bool scaled_imag = pctx.get_scaled(); // may be changed by UCS_string()
const UCS_string ucs_i(value2.cval_i, scaled_imag, pctx);
ucs.append(ucs_i);
info.imag_len = ucs.size() - info.real_len;
if (scaled_imag) info.flags |= imag_has_E;
}
return PrintBuffer(ucs, info);
}
//-----------------------------------------------------------------------------
bool
ComplexCell::need_scaling(const PrintContext &pctx) const
{
// a complex number needs scaling if the real part needs it, ot
// the complex part is significant and needs it.
return FloatCell::need_scaling(value.cval_r, pctx.get_PP()) ||
(!is_near_real() &&
FloatCell::need_scaling(value2.cval_i, pctx.get_PP()));
}
//-----------------------------------------------------------------------------
