1 2/* 3 This file is part of GNU APL, a free implementation of the 4 ISO/IEC Standard 13751, "Programming Language APL, Extended" 5 6 Copyright (C) 2008-2017 Dr. Jürgen Sauermann 7 8 This program is free software: you can redistribute it and/or modify 9 it under the terms of the GNU General Public License as published by 10 the Free Software Foundation, either version 3 of the License, or 11 (at your option) any later version. 12 13 This program is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 GNU General Public License for more details. 17 18 You should have received a copy of the GNU General Public License 19 along with this program. If not, see <http://www.gnu.org/licenses/>. 20*/ 21 22#ifndef __CELL_ICC_DEFINED__ 23#define __CELL_ICC_DEFINED__ 24 25//----------------------------------------------------------------------------- 26inline bool 27Cell::tolerantly_equal(APL_Float A, APL_Float B, double C) 28{ 29 // if the signs of A and B differ then they are unequal (ISO standard 30 // page 19). We treat exact 0.0 as having both signs 31 // 32 if (A == B) return true; 33 if (A < 0.0 && B > 0.0) return false; 34 if (A > 0.0 && B < 0.0) return false; 35 36const APL_Float mag_A = A < 0.0 ? APL_Float(-A) : A; 37const APL_Float mag_B = B < 0.0 ? APL_Float(-B) : B; 38const APL_Float mag_max = mag_A > mag_B ? mag_A : mag_B; 39 40const APL_Float dist_A_B = A > B ? (A - B) : (B - A); 41 42 return dist_A_B < C*mag_max; 43} 44//----------------------------------------------------------------------------- 45inline bool 46Cell::same_half_plane(APL_Complex A, APL_Complex B) 47{ 48 if (A.real() > 0.0 && B.real() > 0.0) return true; 49 if (A.real() < 0.0 && B.real() < 0.0) return true; 50 51 if (A.imag() > 0.0 && B.imag() > 0.0) return true; 52 if (A.imag() < 0.0 && B.imag() < 0.0) return true; 53 54 return false; 55} 56//----------------------------------------------------------------------------- 57inline bool 58Cell::tolerantly_equal(APL_Complex A, APL_Complex B, double C) 59{ 60 // if A equals B, return true 61 // 62 if (A == B) return true; 63 64 // ISO: if A and B are not in the same half-plane, return false. 65 // 66 // Note: this leads to problems with complex modulo of Gaussian integers, so 67 // we do not do it. 68 // 69 // if (!same_half_plane(A, B)) return false; 70 71 // If the distance-between A and B is ≤ C times the larger-magnitude 72 // of A and B, return true 73 // 74 // Implementation: Instead of mag(A-B) ≤ C × max(mag(A), mag(B)) 75 // we compute: mag²(A-B) ≤ C² × max(mag²(A), mag²(B)) 76 // == max(mag²(C×A), mag²(C×B)) 77 // 78 // that avoids unneccessary computations of sqrt(). 79 // 80 81 // 1. compute left side: mag²(A-B) 82 // 83const APL_Complex A_B = A - B; 84const APL_Float dist2_A_B = A_B.real() * A_B.real() 85 + A_B.imag() * A_B.imag(); 86 87 // 2. compute right side: max(mag²C×A, mag²C×B) 88 // 89const APL_Float Cmag2_A = C * ComplexCell::mag2(A); 90const APL_Float Cmag2_B = C * ComplexCell::mag2(B); 91const APL_Float mag2_max = Cmag2_A > Cmag2_B ? Cmag2_A : Cmag2_B; 92 93 // compare 94 // 95 return dist2_A_B <= mag2_max; 96} 97//----------------------------------------------------------------------------- 98inline bool 99Cell::integral_within(APL_Float A, double qct) 100{ 101 if (A < (floor(A) + qct)) return true; 102 if (A > (ceil(A) - qct)) return true; 103 return false; 104 105 // Note: despite its name (coming from the ISO standard) this implementation 106 // of integral_within is NOT the one described in the ISO standard. The 107 // description in the standard was leading to problems in bif_residue 108 // 109 // 1, scale qct to the magnitude of A 110 // 111 qct = A * A * qct; 112 113const APL_Float real_up = ceil (A) - qct; 114const APL_Float real_dn = floor(A) + qct; 115 116 return A > real_up || A < real_dn; 117} 118//----------------------------------------------------------------------------- 119 120#endif // __CELL_ICC_DEFINED__ 121 122