From cdb84e8f913b5573da8ccef12b88ac48c09cf54e Mon Sep 17 00:00:00 2001 From: baudin Date: Fri, 4 Jun 2010 12:56:48 +0000 Subject: Fixed missing copyright notice. --- src/c/elementaryFunctions/acos/cacoss.c | 293 ++++++++-------- src/c/elementaryFunctions/acos/zacoss.c | 293 ++++++++-------- src/c/elementaryFunctions/asin/casins.c | 290 ++++++++-------- src/c/elementaryFunctions/asin/zasins.c | 290 ++++++++-------- src/c/elementaryFunctions/atan/catans.c | 497 ++++++++++++++-------------- src/c/elementaryFunctions/atan/zatans.c | 483 +++++++++++++-------------- src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c | 153 ++++----- src/c/elementaryFunctions/lnp1m1/slnp1m1s.c | 153 ++++----- src/c/elementaryFunctions/log/clogs.c | 129 ++++---- src/c/elementaryFunctions/log/zlogs.c | 129 ++++---- src/c/elementaryFunctions/log1p/dlog1ps.c | 67 ++-- src/c/elementaryFunctions/log1p/slog1ps.c | 68 ++-- src/c/elementaryFunctions/sqrt/csqrts.c | 1 + src/c/elementaryFunctions/sqrt/zsqrts.c | 1 + src/c/elementaryFunctions/tan/ztans.c | 207 ++++++------ 15 files changed, 1536 insertions(+), 1518 deletions(-) (limited to 'src') diff --git a/src/c/elementaryFunctions/acos/cacoss.c b/src/c/elementaryFunctions/acos/cacoss.c index 1059c8c3..97420313 100644 --- a/src/c/elementaryFunctions/acos/cacoss.c +++ b/src/c/elementaryFunctions/acos/cacoss.c @@ -1,146 +1,147 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -/* - * This fonction is a translation of fortran wacos write by Bruno Pincon - * REFERENCE - * This is a Fortran-77 translation of an algorithm by - * T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which - * appears in their article : - * "Implementing the Complex Arcsine and Arccosine - * Functions Using Exception Handling", ACM, TOMS, - * Vol 23, No. 3, Sept 1997, p. 299-335 - */ - -#include "acos.h" -#include "atan.h" -#include "log.h" -#include "log1p.h" -#include "sqrt.h" -#include "abs.h" -#include "lapack.h" -#include "min.h" -#include "max.h" - -#define localSign(x) (x>0 ? 1.0f : -1.0f) - -floatComplex cacoss(floatComplex z) { - static float sfltPi = 3.1415926535897932384626433f; - static float sfltPi_2 = 1.5707963267948966192313216f; - static float sfltLn2 = 0.6931471805599453094172321f; - static float sfltAcross = 1.5f; - static float sfltBcross = 0.6417f; - - float fltLsup = ssqrts((float) getOverflowThreshold())/8.0f; - float fltLinf = 4.0f * ssqrts((float) getUnderflowThreshold()); - float fltEpsm = ssqrts((float) getRelativeMachinePrecision()); - - float fltAbsReal = sabss(creals(z)); - float fltAbsImg = sabss(cimags(z)); - float fltSignReal = localSign(creals(z)); - float fltSignImg = localSign(cimags(z)); - - float fltR = 0, fltS = 0, fltA = 0, fltB = 0; - - float fltTemp = 0; - - float _pfltReal = 0; - float _pfltImg = 0; - - if( min(fltAbsReal, fltAbsImg) > fltLinf && max(fltAbsReal, fltAbsImg) <= fltLsup) - {/* we are in the safe region */ - fltR = ssqrts( (fltAbsReal + 1 )*(fltAbsReal + 1 ) + fltAbsImg*fltAbsImg); - fltS = ssqrts( (fltAbsReal - 1 )*(fltAbsReal - 1 ) + fltAbsImg*fltAbsImg); - fltA = 0.5f * ( fltR + fltS ); - fltB = fltAbsReal / fltA; - - - /* compute the real part */ - if(fltB <= sfltBcross) - _pfltReal = sacoss(fltB); - else if( fltAbsReal <= 1) - _pfltReal = satans(ssqrts(0.5f * (fltA + fltAbsReal) * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + (fltS + (1 - fltAbsReal)))) / fltAbsReal); - else - _pfltReal = satans((fltAbsImg * ssqrts(0.5f * ((fltA + fltAbsReal) / (fltR + (fltAbsReal + 1)) + (fltA + fltAbsReal) / (fltS + (fltAbsReal - 1))))) / fltAbsReal); - - /* compute the imaginary part */ - if(fltA <= sfltAcross) - { - float fltImg1 = 0; - - if(fltAbsReal < 1) - /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */ - fltImg1 = 0.5f * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + fltAbsImg*fltAbsImg / (fltS + (1 - fltAbsReal))); - else - /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */ - fltImg1 = 0.5f * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + (fltS + (fltAbsReal - 1))); - /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */ - fltTemp = fltImg1 + ssqrts(fltImg1 *( fltA + 1)); - _pfltImg = slog1ps(fltTemp); - } - else - /* ai = log(A + sqrt(A**2 - 1.d0)) */ - _pfltImg = slogs(fltA + ssqrts(fltA*fltA - 1)); - } - else - {/* evaluation in the special regions ... */ - if(fltAbsImg <= fltEpsm * sabss(fltAbsReal - 1)) - { - if(fltAbsReal < 1) - { - _pfltReal = sacoss(fltAbsReal); - _pfltImg = fltAbsImg / ssqrts((1 + fltAbsReal) * (1 - fltAbsReal)); - } - else - { - _pfltReal = 0; - if(fltAbsReal <= fltLsup) - { - fltTemp = (fltAbsReal - 1) + ssqrts((fltAbsReal - 1) * (fltAbsReal + 1)); - _pfltImg = slog1ps(fltTemp); - } - else - _pfltImg = sfltLn2 + slogs(fltAbsReal); - } - } - else if(fltAbsImg < fltLinf) - { - _pfltReal = ssqrts(fltAbsImg); - _pfltImg = _pfltReal; - } - else if((fltEpsm * fltAbsImg - 1 >= fltAbsReal)) - { - _pfltReal = sfltPi_2; - _pfltImg = sfltLn2 + slogs(fltAbsImg); - } - else if(fltAbsReal > 1) - { - _pfltReal = satans(fltAbsImg / fltAbsReal); - fltTemp = (fltAbsReal / fltAbsImg)*(fltAbsReal / fltAbsImg); - _pfltImg = sfltLn2 + slogs(fltAbsImg) + 0.5f * slog1ps(fltTemp); - } - else - { - float fltTemp2 = ssqrts(1 + fltAbsImg*fltAbsImg); - _pfltReal = sfltPi_2; - fltTemp = 2 * fltAbsImg * (fltAbsImg + fltTemp2); - _pfltImg = 0.5f * slog1ps(fltTemp); - } - } - if(fltSignReal < 0) - _pfltReal = sfltPi - _pfltReal; - - if(fltAbsImg != 0 || fltSignReal < 0) - _pfltImg = - fltSignImg * _pfltImg; - - return FloatComplex(_pfltReal, _pfltImg); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +/* + * This fonction is a translation of fortran wacos write by Bruno Pincon + * REFERENCE + * This is a Fortran-77 translation of an algorithm by + * T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which + * appears in their article : + * "Implementing the Complex Arcsine and Arccosine + * Functions Using Exception Handling", ACM, TOMS, + * Vol 23, No. 3, Sept 1997, p. 299-335 + */ + +#include "acos.h" +#include "atan.h" +#include "log.h" +#include "log1p.h" +#include "sqrt.h" +#include "abs.h" +#include "lapack.h" +#include "min.h" +#include "max.h" + +#define localSign(x) (x>0 ? 1.0f : -1.0f) + +floatComplex cacoss(floatComplex z) { + static float sfltPi = 3.1415926535897932384626433f; + static float sfltPi_2 = 1.5707963267948966192313216f; + static float sfltLn2 = 0.6931471805599453094172321f; + static float sfltAcross = 1.5f; + static float sfltBcross = 0.6417f; + + float fltLsup = ssqrts((float) getOverflowThreshold())/8.0f; + float fltLinf = 4.0f * ssqrts((float) getUnderflowThreshold()); + float fltEpsm = ssqrts((float) getRelativeMachinePrecision()); + + float fltAbsReal = sabss(creals(z)); + float fltAbsImg = sabss(cimags(z)); + float fltSignReal = localSign(creals(z)); + float fltSignImg = localSign(cimags(z)); + + float fltR = 0, fltS = 0, fltA = 0, fltB = 0; + + float fltTemp = 0; + + float _pfltReal = 0; + float _pfltImg = 0; + + if( min(fltAbsReal, fltAbsImg) > fltLinf && max(fltAbsReal, fltAbsImg) <= fltLsup) + {/* we are in the safe region */ + fltR = ssqrts( (fltAbsReal + 1 )*(fltAbsReal + 1 ) + fltAbsImg*fltAbsImg); + fltS = ssqrts( (fltAbsReal - 1 )*(fltAbsReal - 1 ) + fltAbsImg*fltAbsImg); + fltA = 0.5f * ( fltR + fltS ); + fltB = fltAbsReal / fltA; + + + /* compute the real part */ + if(fltB <= sfltBcross) + _pfltReal = sacoss(fltB); + else if( fltAbsReal <= 1) + _pfltReal = satans(ssqrts(0.5f * (fltA + fltAbsReal) * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + (fltS + (1 - fltAbsReal)))) / fltAbsReal); + else + _pfltReal = satans((fltAbsImg * ssqrts(0.5f * ((fltA + fltAbsReal) / (fltR + (fltAbsReal + 1)) + (fltA + fltAbsReal) / (fltS + (fltAbsReal - 1))))) / fltAbsReal); + + /* compute the imaginary part */ + if(fltA <= sfltAcross) + { + float fltImg1 = 0; + + if(fltAbsReal < 1) + /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */ + fltImg1 = 0.5f * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + fltAbsImg*fltAbsImg / (fltS + (1 - fltAbsReal))); + else + /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */ + fltImg1 = 0.5f * (fltAbsImg*fltAbsImg / (fltR + (fltAbsReal + 1)) + (fltS + (fltAbsReal - 1))); + /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */ + fltTemp = fltImg1 + ssqrts(fltImg1 *( fltA + 1)); + _pfltImg = slog1ps(fltTemp); + } + else + /* ai = log(A + sqrt(A**2 - 1.d0)) */ + _pfltImg = slogs(fltA + ssqrts(fltA*fltA - 1)); + } + else + {/* evaluation in the special regions ... */ + if(fltAbsImg <= fltEpsm * sabss(fltAbsReal - 1)) + { + if(fltAbsReal < 1) + { + _pfltReal = sacoss(fltAbsReal); + _pfltImg = fltAbsImg / ssqrts((1 + fltAbsReal) * (1 - fltAbsReal)); + } + else + { + _pfltReal = 0; + if(fltAbsReal <= fltLsup) + { + fltTemp = (fltAbsReal - 1) + ssqrts((fltAbsReal - 1) * (fltAbsReal + 1)); + _pfltImg = slog1ps(fltTemp); + } + else + _pfltImg = sfltLn2 + slogs(fltAbsReal); + } + } + else if(fltAbsImg < fltLinf) + { + _pfltReal = ssqrts(fltAbsImg); + _pfltImg = _pfltReal; + } + else if((fltEpsm * fltAbsImg - 1 >= fltAbsReal)) + { + _pfltReal = sfltPi_2; + _pfltImg = sfltLn2 + slogs(fltAbsImg); + } + else if(fltAbsReal > 1) + { + _pfltReal = satans(fltAbsImg / fltAbsReal); + fltTemp = (fltAbsReal / fltAbsImg)*(fltAbsReal / fltAbsImg); + _pfltImg = sfltLn2 + slogs(fltAbsImg) + 0.5f * slog1ps(fltTemp); + } + else + { + float fltTemp2 = ssqrts(1 + fltAbsImg*fltAbsImg); + _pfltReal = sfltPi_2; + fltTemp = 2 * fltAbsImg * (fltAbsImg + fltTemp2); + _pfltImg = 0.5f * slog1ps(fltTemp); + } + } + if(fltSignReal < 0) + _pfltReal = sfltPi - _pfltReal; + + if(fltAbsImg != 0 || fltSignReal < 0) + _pfltImg = - fltSignImg * _pfltImg; + + return FloatComplex(_pfltReal, _pfltImg); +} diff --git a/src/c/elementaryFunctions/acos/zacoss.c b/src/c/elementaryFunctions/acos/zacoss.c index 7758a932..de6f3fe9 100644 --- a/src/c/elementaryFunctions/acos/zacoss.c +++ b/src/c/elementaryFunctions/acos/zacoss.c @@ -1,146 +1,147 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -/* - * This fonction is a translation of fortran wacos write by Bruno Pincon - * REFERENCE - * This is a Fortran-77 translation of an algorithm by - * T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which - * appears in their article : - * "Implementing the Complex Arcsine and Arccosine - * Functions Using Exception Handling", ACM, TOMS, - * Vol 23, No. 3, Sept 1997, p. 299-335 - */ - -#include "acos.h" -#include "atan.h" -#include "log.h" -#include "log1p.h" -#include "sqrt.h" -#include "abs.h" -#include "lapack.h" -#include "min.h" -#include "max.h" - -#define localSign(x) (x>0 ? 1 : -1) - -doubleComplex zacoss(doubleComplex z) { - static double sdblPi = 3.1415926535897932384626433; - static double sdblPi_2 = 1.5707963267948966192313216; - static double sdblLn2 = 0.6931471805599453094172321; - static double sdblAcross = 1.5; - static double sdblBcross = 0.6417; - - double dblLsup = dsqrts(getOverflowThreshold())/8.0; - double dblLinf = 4.0 * dsqrts(getUnderflowThreshold()); - double dblEpsm = dsqrts(getRelativeMachinePrecision()); - - double dblAbsReal = dabss(zreals(z)); - double dblAbsImg = dabss(zimags(z)); - double dblSignReal = localSign(zreals(z)); - double dblSignImg = localSign(zimags(z)); - - double dblR = 0, dblS = 0, dblA = 0, dblB = 0; - - double dblTemp = 0; - - double _pdblReal = 0; - double _pdblImg = 0; - - if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup) - {/* we are in the safe region */ - dblR = dsqrts( (dblAbsReal + 1 )*(dblAbsReal + 1 ) + dblAbsImg*dblAbsImg); - dblS = dsqrts( (dblAbsReal - 1 )*(dblAbsReal - 1 ) + dblAbsImg*dblAbsImg); - dblA = 0.5 * ( dblR + dblS ); - dblB = dblAbsReal / dblA; - - - /* compute the real part */ - if(dblB <= sdblBcross) - _pdblReal = dacoss(dblB); - else if( dblAbsReal <= 1) - _pdblReal = datans(dsqrts(0.5 * (dblA + dblAbsReal) * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal)))) / dblAbsReal); - else - _pdblReal = datans((dblAbsImg * dsqrts(0.5 * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal - 1))))) / dblAbsReal); - - /* compute the imaginary part */ - if(dblA <= sdblAcross) - { - double dblImg1 = 0; - - if(dblAbsReal < 1) - /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */ - dblImg1 = 0.5 * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg*dblAbsImg / (dblS + (1 - dblAbsReal))); - else - /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */ - dblImg1 = 0.5 * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1))); - /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */ - dblTemp = dblImg1 + dsqrts(dblImg1 *( dblA + 1)); - _pdblImg = dlog1ps(dblTemp); - } - else - /* ai = log(A + sqrt(A**2 - 1.d0)) */ - _pdblImg = dlogs(dblA + dsqrts(dblA*dblA - 1)); - } - else - {/* evaluation in the special regions ... */ - if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1)) - { - if(dblAbsReal < 1) - { - _pdblReal = dacoss(dblAbsReal); - _pdblImg = dblAbsImg / dsqrts((1 + dblAbsReal) * (1 - dblAbsReal)); - } - else - { - _pdblReal = 0; - if(dblAbsReal <= dblLsup) - { - dblTemp = (dblAbsReal - 1) + dsqrts((dblAbsReal - 1) * (dblAbsReal + 1)); - _pdblImg = dlog1ps(dblTemp); - } - else - _pdblImg = sdblLn2 + dlogs(dblAbsReal); - } - } - else if(dblAbsImg < dblLinf) - { - _pdblReal = dsqrts(dblAbsImg); - _pdblImg = _pdblReal; - } - else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal)) - { - _pdblReal = sdblPi_2; - _pdblImg = sdblLn2 + dlogs(dblAbsImg); - } - else if(dblAbsReal > 1) - { - _pdblReal = datans(dblAbsImg / dblAbsReal); - dblTemp = (dblAbsReal / dblAbsImg)*(dblAbsReal / dblAbsImg); - _pdblImg = sdblLn2 + dlogs(dblAbsImg) + 0.5 * dlog1ps(dblTemp); - } - else - { - double dblTemp2 = dsqrts(1 + dblAbsImg*dblAbsImg); - _pdblReal = sdblPi_2; - dblTemp = 2 * dblAbsImg * (dblAbsImg + dblTemp2); - _pdblImg = 0.5 * dlog1ps(dblTemp); - } - } - if(dblSignReal < 0) - _pdblReal = sdblPi - _pdblReal; - - if(dblAbsImg != 0 || dblSignReal < 0) - _pdblImg = - dblSignImg * _pdblImg; - - return DoubleComplex(_pdblReal, _pdblImg); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +/* + * This fonction is a translation of fortran wacos write by Bruno Pincon + * REFERENCE + * This is a Fortran-77 translation of an algorithm by + * T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which + * appears in their article : + * "Implementing the Complex Arcsine and Arccosine + * Functions Using Exception Handling", ACM, TOMS, + * Vol 23, No. 3, Sept 1997, p. 299-335 + */ + +#include "acos.h" +#include "atan.h" +#include "log.h" +#include "log1p.h" +#include "sqrt.h" +#include "abs.h" +#include "lapack.h" +#include "min.h" +#include "max.h" + +#define localSign(x) (x>0 ? 1 : -1) + +doubleComplex zacoss(doubleComplex z) { + static double sdblPi = 3.1415926535897932384626433; + static double sdblPi_2 = 1.5707963267948966192313216; + static double sdblLn2 = 0.6931471805599453094172321; + static double sdblAcross = 1.5; + static double sdblBcross = 0.6417; + + double dblLsup = dsqrts(getOverflowThreshold())/8.0; + double dblLinf = 4.0 * dsqrts(getUnderflowThreshold()); + double dblEpsm = dsqrts(getRelativeMachinePrecision()); + + double dblAbsReal = dabss(zreals(z)); + double dblAbsImg = dabss(zimags(z)); + double dblSignReal = localSign(zreals(z)); + double dblSignImg = localSign(zimags(z)); + + double dblR = 0, dblS = 0, dblA = 0, dblB = 0; + + double dblTemp = 0; + + double _pdblReal = 0; + double _pdblImg = 0; + + if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup) + {/* we are in the safe region */ + dblR = dsqrts( (dblAbsReal + 1 )*(dblAbsReal + 1 ) + dblAbsImg*dblAbsImg); + dblS = dsqrts( (dblAbsReal - 1 )*(dblAbsReal - 1 ) + dblAbsImg*dblAbsImg); + dblA = 0.5 * ( dblR + dblS ); + dblB = dblAbsReal / dblA; + + + /* compute the real part */ + if(dblB <= sdblBcross) + _pdblReal = dacoss(dblB); + else if( dblAbsReal <= 1) + _pdblReal = datans(dsqrts(0.5 * (dblA + dblAbsReal) * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal)))) / dblAbsReal); + else + _pdblReal = datans((dblAbsImg * dsqrts(0.5 * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal - 1))))) / dblAbsReal); + + /* compute the imaginary part */ + if(dblA <= sdblAcross) + { + double dblImg1 = 0; + + if(dblAbsReal < 1) + /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */ + dblImg1 = 0.5 * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg*dblAbsImg / (dblS + (1 - dblAbsReal))); + else + /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */ + dblImg1 = 0.5 * (dblAbsImg*dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1))); + /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */ + dblTemp = dblImg1 + dsqrts(dblImg1 *( dblA + 1)); + _pdblImg = dlog1ps(dblTemp); + } + else + /* ai = log(A + sqrt(A**2 - 1.d0)) */ + _pdblImg = dlogs(dblA + dsqrts(dblA*dblA - 1)); + } + else + {/* evaluation in the special regions ... */ + if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1)) + { + if(dblAbsReal < 1) + { + _pdblReal = dacoss(dblAbsReal); + _pdblImg = dblAbsImg / dsqrts((1 + dblAbsReal) * (1 - dblAbsReal)); + } + else + { + _pdblReal = 0; + if(dblAbsReal <= dblLsup) + { + dblTemp = (dblAbsReal - 1) + dsqrts((dblAbsReal - 1) * (dblAbsReal + 1)); + _pdblImg = dlog1ps(dblTemp); + } + else + _pdblImg = sdblLn2 + dlogs(dblAbsReal); + } + } + else if(dblAbsImg < dblLinf) + { + _pdblReal = dsqrts(dblAbsImg); + _pdblImg = _pdblReal; + } + else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal)) + { + _pdblReal = sdblPi_2; + _pdblImg = sdblLn2 + dlogs(dblAbsImg); + } + else if(dblAbsReal > 1) + { + _pdblReal = datans(dblAbsImg / dblAbsReal); + dblTemp = (dblAbsReal / dblAbsImg)*(dblAbsReal / dblAbsImg); + _pdblImg = sdblLn2 + dlogs(dblAbsImg) + 0.5 * dlog1ps(dblTemp); + } + else + { + double dblTemp2 = dsqrts(1 + dblAbsImg*dblAbsImg); + _pdblReal = sdblPi_2; + dblTemp = 2 * dblAbsImg * (dblAbsImg + dblTemp2); + _pdblImg = 0.5 * dlog1ps(dblTemp); + } + } + if(dblSignReal < 0) + _pdblReal = sdblPi - _pdblReal; + + if(dblAbsImg != 0 || dblSignReal < 0) + _pdblImg = - dblSignImg * _pdblImg; + + return DoubleComplex(_pdblReal, _pdblImg); +} diff --git a/src/c/elementaryFunctions/asin/casins.c b/src/c/elementaryFunctions/asin/casins.c index 9c15ce53..35a4a8d8 100644 --- a/src/c/elementaryFunctions/asin/casins.c +++ b/src/c/elementaryFunctions/asin/casins.c @@ -1,144 +1,146 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -/* - * REFERENCE - * This is a Fortran-77 translation of an algorithm by - * T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which - * appears in their article : - * "Implementing the Complex Arcsine and Arccosine - * Functions Using Exception Handling", ACM, TOMS, - * Vol 23, No. 3, Sept 1997, p. 299-335 - */ - -#include "lapack.h" -#include "asin.h" -#include "atan.h" -#include "sqrt.h" -#include "abs.h" -#include "log.h" -#include "log1p.h" -#include "min.h" -#include "max.h" - -floatComplex casins(floatComplex z) { - static float sdblPi_2 = 1.5707963267948966192313216f; - static float sdblLn2 = 0.6931471805599453094172321f; - static float sdblAcross = 1.5f; - static float sdblBcross = 0.6417f; - - float dblLsup = ssqrts((float) getOverflowThreshold())/ 8.0f; - float dblLinf = 4.0f * ssqrts((float) getUnderflowThreshold()); - float dblEpsm = ssqrts((float) getRelativeMachinePrecision()); - - float _dblReal = creals(z); - float _dblImg = cimags(z); - - float dblAbsReal = sabss(_dblReal); - float dblAbsImg = sabss(_dblImg); - float iSignReal = _dblReal < 0 ? -1.0f : 1.0f; - float iSignImg = _dblImg < 0 ? -1.0f : 1.0f; - - float dblR = 0, dblS = 0, dblA = 0, dblB = 0; - - float dblTemp = 0; - - float _pdblReal = 0; - float _pdblImg = 0; - - if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup) - { - /* we are in the safe region */ - dblR = ssqrts( (dblAbsReal + 1) * (dblAbsReal + 1) + dblAbsImg * dblAbsImg); - dblS = ssqrts( (dblAbsReal - 1) * (dblAbsReal - 1) + dblAbsImg * dblAbsImg); - dblA = (float) 0.5 * ( dblR + dblS ); - dblB = dblAbsReal / dblA; - - - /* compute the real part */ - if(dblB <= sdblBcross) - _pdblReal = sasins(dblB); - else if( dblAbsReal <= 1) - _pdblReal = satans(dblAbsReal / ssqrts( 0.5f * (dblA + dblAbsReal) * ( (dblAbsImg * dblAbsImg) / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal))))); - else - _pdblReal = satans(dblAbsReal / (dblAbsImg * ssqrts( 0.5f * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal-1)))))); - - /* compute the imaginary part */ - if(dblA <= sdblAcross) - { - float dblImg1 = 0; - - if(dblAbsReal < 1) - /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */ - dblImg1 = 0.5f * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg * dblAbsImg / (dblS + (1 - dblAbsReal))); - else - /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */ - dblImg1 = 0.5f * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1))); - /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */ - dblTemp = dblImg1 + ssqrts(dblImg1 * (dblA + 1)); - _pdblImg = slog1ps(dblTemp); - } - else - /* ai = log(A + sqrt(A**2 - 1.d0)) */ - _pdblImg = slogs(dblA + ssqrts(dblA * dblA - (float) 1.0)); - } - else - { - /* evaluation in the special regions ... */ - if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1)) - { - if(dblAbsReal < 1) - { - _pdblReal = sasins(dblAbsReal); - _pdblImg = dblAbsImg / ssqrts((1 + dblAbsReal) * (1 - dblAbsReal)); - } - else - { - _pdblReal = sdblPi_2; - if(dblAbsReal <= dblLsup) - { - dblTemp = (dblAbsReal - 1) + ssqrts((dblAbsReal - 1) * (dblAbsReal + 1)); - _pdblImg = slog1ps(dblTemp); - } - else - _pdblImg = sdblLn2 + slogs(dblAbsReal); - } - } - else if(dblAbsImg < dblLinf) - { - _pdblReal = sdblPi_2 - ssqrts(dblAbsImg); - _pdblImg = ssqrts(dblAbsImg); - } - else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal)) - { - _pdblReal = dblAbsReal * dblAbsImg; - _pdblImg = sdblLn2 + slogs(dblAbsReal); - } - else if(dblAbsReal > 1) - { - _pdblReal = satans(dblAbsReal / dblAbsImg); - dblTemp = (dblAbsReal / dblAbsImg) * (dblAbsReal / dblAbsImg); - _pdblImg = sdblLn2 + slogs(dblAbsReal) + 0.5f * slog1ps(dblTemp); - } - else - { - float dblTemp2 = ssqrts(1 + dblAbsImg * dblAbsImg); - _pdblReal = dblAbsReal / dblTemp2; - dblTemp = 2.0f * dblAbsImg * (dblAbsImg + dblTemp2); - _pdblImg = 0.5f * slog1ps(dblTemp); - } - } - _pdblReal *= iSignReal; - _pdblImg *= iSignImg; - - return (FloatComplex(_pdblReal, _pdblImg)); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +/* + * REFERENCE + * This is a Fortran-77 translation of an algorithm by + * T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which + * appears in their article : + * "Implementing the Complex Arcsine and Arccosine + * Functions Using Exception Handling", ACM, TOMS, + * Vol 23, No. 3, Sept 1997, p. 299-335 + * Thanks to Tom Fairgrieve + */ + +#include "lapack.h" +#include "asin.h" +#include "atan.h" +#include "sqrt.h" +#include "abs.h" +#include "log.h" +#include "log1p.h" +#include "min.h" +#include "max.h" + +floatComplex casins(floatComplex z) { + static float sdblPi_2 = 1.5707963267948966192313216f; + static float sdblLn2 = 0.6931471805599453094172321f; + static float sdblAcross = 1.5f; + static float sdblBcross = 0.6417f; + + float dblLsup = ssqrts((float) getOverflowThreshold())/ 8.0f; + float dblLinf = 4.0f * ssqrts((float) getUnderflowThreshold()); + float dblEpsm = ssqrts((float) getRelativeMachinePrecision()); + + float _dblReal = creals(z); + float _dblImg = cimags(z); + + float dblAbsReal = sabss(_dblReal); + float dblAbsImg = sabss(_dblImg); + float iSignReal = _dblReal < 0 ? -1.0f : 1.0f; + float iSignImg = _dblImg < 0 ? -1.0f : 1.0f; + + float dblR = 0, dblS = 0, dblA = 0, dblB = 0; + + float dblTemp = 0; + + float _pdblReal = 0; + float _pdblImg = 0; + + if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup) + { + /* we are in the safe region */ + dblR = ssqrts( (dblAbsReal + 1) * (dblAbsReal + 1) + dblAbsImg * dblAbsImg); + dblS = ssqrts( (dblAbsReal - 1) * (dblAbsReal - 1) + dblAbsImg * dblAbsImg); + dblA = (float) 0.5 * ( dblR + dblS ); + dblB = dblAbsReal / dblA; + + + /* compute the real part */ + if(dblB <= sdblBcross) + _pdblReal = sasins(dblB); + else if( dblAbsReal <= 1) + _pdblReal = satans(dblAbsReal / ssqrts( 0.5f * (dblA + dblAbsReal) * ( (dblAbsImg * dblAbsImg) / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal))))); + else + _pdblReal = satans(dblAbsReal / (dblAbsImg * ssqrts( 0.5f * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal-1)))))); + + /* compute the imaginary part */ + if(dblA <= sdblAcross) + { + float dblImg1 = 0; + + if(dblAbsReal < 1) + /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */ + dblImg1 = 0.5f * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg * dblAbsImg / (dblS + (1 - dblAbsReal))); + else + /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */ + dblImg1 = 0.5f * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1))); + /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */ + dblTemp = dblImg1 + ssqrts(dblImg1 * (dblA + 1)); + _pdblImg = slog1ps(dblTemp); + } + else + /* ai = log(A + sqrt(A**2 - 1.d0)) */ + _pdblImg = slogs(dblA + ssqrts(dblA * dblA - (float) 1.0)); + } + else + { + /* evaluation in the special regions ... */ + if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1)) + { + if(dblAbsReal < 1) + { + _pdblReal = sasins(dblAbsReal); + _pdblImg = dblAbsImg / ssqrts((1 + dblAbsReal) * (1 - dblAbsReal)); + } + else + { + _pdblReal = sdblPi_2; + if(dblAbsReal <= dblLsup) + { + dblTemp = (dblAbsReal - 1) + ssqrts((dblAbsReal - 1) * (dblAbsReal + 1)); + _pdblImg = slog1ps(dblTemp); + } + else + _pdblImg = sdblLn2 + slogs(dblAbsReal); + } + } + else if(dblAbsImg < dblLinf) + { + _pdblReal = sdblPi_2 - ssqrts(dblAbsImg); + _pdblImg = ssqrts(dblAbsImg); + } + else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal)) + { + _pdblReal = dblAbsReal * dblAbsImg; + _pdblImg = sdblLn2 + slogs(dblAbsReal); + } + else if(dblAbsReal > 1) + { + _pdblReal = satans(dblAbsReal / dblAbsImg); + dblTemp = (dblAbsReal / dblAbsImg) * (dblAbsReal / dblAbsImg); + _pdblImg = sdblLn2 + slogs(dblAbsReal) + 0.5f * slog1ps(dblTemp); + } + else + { + float dblTemp2 = ssqrts(1 + dblAbsImg * dblAbsImg); + _pdblReal = dblAbsReal / dblTemp2; + dblTemp = 2.0f * dblAbsImg * (dblAbsImg + dblTemp2); + _pdblImg = 0.5f * slog1ps(dblTemp); + } + } + _pdblReal *= iSignReal; + _pdblImg *= iSignImg; + + return (FloatComplex(_pdblReal, _pdblImg)); +} diff --git a/src/c/elementaryFunctions/asin/zasins.c b/src/c/elementaryFunctions/asin/zasins.c index 5d2110ff..5bd586a8 100644 --- a/src/c/elementaryFunctions/asin/zasins.c +++ b/src/c/elementaryFunctions/asin/zasins.c @@ -1,144 +1,146 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -/* - * REFERENCE - * This is a Fortran-77 translation of an algorithm by - * T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which - * appears in their article : - * "Implementing the Complex Arcsine and Arccosine - * Functions Using Exception Handling", ACM, TOMS, - * Vol 23, No. 3, Sept 1997, p. 299-335 - */ - -#include "lapack.h" -#include "asin.h" -#include "atan.h" -#include "sqrt.h" -#include "abs.h" -#include "log.h" -#include "log1p.h" -#include "min.h" -#include "max.h" - -doubleComplex zasins(doubleComplex z) { - static double sdblPi_2 = 1.5707963267948966192313216; - static double sdblLn2 = 0.6931471805599453094172321; - static double sdblAcross = 1.5; - static double sdblBcross = 0.6417; - - double dblLsup = dsqrts(getOverflowThreshold())/8.0; - double dblLinf = 4 * dsqrts(getUnderflowThreshold()); - double dblEpsm = dsqrts(getRelativeMachinePrecision()); - - double _dblReal = zreals(z); - double _dblImg = zimags(z); - - double dblAbsReal = dabss(_dblReal); - double dblAbsImg = dabss(_dblImg); - int iSignReal = _dblReal < 0 ? -1 : 1; - int iSignImg = _dblImg < 0 ? -1 : 1; - - double dblR = 0, dblS = 0, dblA = 0, dblB = 0; - - double dblTemp = 0; - - double _pdblReal = 0; - double _pdblImg = 0; - - if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup) - { - /* we are in the safe region */ - dblR = dsqrts( (dblAbsReal + 1) * (dblAbsReal + 1) + dblAbsImg * dblAbsImg); - dblS = dsqrts( (dblAbsReal - 1) * (dblAbsReal - 1) + dblAbsImg * dblAbsImg); - dblA = 0.5 * ( dblR + dblS ); - dblB = dblAbsReal / dblA; - - - /* compute the real part */ - if(dblB <= sdblBcross) - _pdblReal = dasins(dblB); - else if( dblAbsReal <= 1) - _pdblReal = datans(dblAbsReal / dsqrts( 0.5 * (dblA + dblAbsReal) * ( (dblAbsImg * dblAbsImg) / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal))))); - else - _pdblReal = datans(dblAbsReal / (dblAbsImg * dsqrts(0.5 * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal-1)))))); - - /* compute the imaginary part */ - if(dblA <= sdblAcross) - { - double dblImg1 = 0; - - if(dblAbsReal < 1) - /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */ - dblImg1 = 0.5 * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg * dblAbsImg / (dblS + (1 - dblAbsReal))); - else - /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */ - dblImg1 = 0.5 * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1))); - /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */ - dblTemp = dblImg1 + dsqrts(dblImg1 * (dblA + 1)); - _pdblImg = dlog1ps(dblTemp); - } - else - /* ai = log(A + sqrt(A**2 - 1.d0)) */ - _pdblImg = dlogs(dblA + dsqrts(dblA * dblA - 1)); - } - else - { - /* evaluation in the special regions ... */ - if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1)) - { - if(dblAbsReal < 1) - { - _pdblReal = dasins(dblAbsReal); - _pdblImg = dblAbsImg / dsqrts((1 + dblAbsReal) * (1 - dblAbsReal)); - } - else - { - _pdblReal = sdblPi_2; - if(dblAbsReal <= dblLsup) - { - dblTemp = (dblAbsReal - 1) + dsqrts((dblAbsReal - 1) * (dblAbsReal + 1)); - _pdblImg = dlog1ps(dblTemp); - } - else - _pdblImg = sdblLn2 + dlogs(dblAbsReal); - } - } - else if(dblAbsImg < dblLinf) - { - _pdblReal = sdblPi_2 - dsqrts(dblAbsImg); - _pdblImg = dsqrts(dblAbsImg); - } - else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal)) - { - _pdblReal = dblAbsReal * dblAbsImg; - _pdblImg = sdblLn2 + dlogs(dblAbsReal); - } - else if(dblAbsReal > 1) - { - _pdblReal = datans(dblAbsReal / dblAbsImg); - dblTemp = (dblAbsReal / dblAbsImg) * (dblAbsReal / dblAbsImg); - _pdblImg = sdblLn2 + dlogs(dblAbsReal) + 0.5 * dlog1ps(dblTemp); - } - else - { - double dblTemp2 = dsqrts(1 + dblAbsImg * dblAbsImg); - _pdblReal = dblAbsReal / dblTemp2; - dblTemp = 2 * dblAbsImg * (dblAbsImg + dblTemp2); - _pdblImg = 0.5 * dlog1ps(dblTemp); - } - } - _pdblReal *= iSignReal; - _pdblImg *= iSignImg; - - return (DoubleComplex(_pdblReal, _pdblImg)); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +/* + * REFERENCE + * This is a Fortran-77 translation of an algorithm by + * T.E. Hull, T. F. Fairgrieve and P.T.P. Tang which + * appears in their article : + * "Implementing the Complex Arcsine and Arccosine + * Functions Using Exception Handling", ACM, TOMS, + * Vol 23, No. 3, Sept 1997, p. 299-335 + * Thanks to Tom Fairgrieve + */ + +#include "lapack.h" +#include "asin.h" +#include "atan.h" +#include "sqrt.h" +#include "abs.h" +#include "log.h" +#include "log1p.h" +#include "min.h" +#include "max.h" + +doubleComplex zasins(doubleComplex z) { + static double sdblPi_2 = 1.5707963267948966192313216; + static double sdblLn2 = 0.6931471805599453094172321; + static double sdblAcross = 1.5; + static double sdblBcross = 0.6417; + + double dblLsup = dsqrts(getOverflowThreshold())/8.0; + double dblLinf = 4 * dsqrts(getUnderflowThreshold()); + double dblEpsm = dsqrts(getRelativeMachinePrecision()); + + double _dblReal = zreals(z); + double _dblImg = zimags(z); + + double dblAbsReal = dabss(_dblReal); + double dblAbsImg = dabss(_dblImg); + int iSignReal = _dblReal < 0 ? -1 : 1; + int iSignImg = _dblImg < 0 ? -1 : 1; + + double dblR = 0, dblS = 0, dblA = 0, dblB = 0; + + double dblTemp = 0; + + double _pdblReal = 0; + double _pdblImg = 0; + + if( min(dblAbsReal, dblAbsImg) > dblLinf && max(dblAbsReal, dblAbsImg) <= dblLsup) + { + /* we are in the safe region */ + dblR = dsqrts( (dblAbsReal + 1) * (dblAbsReal + 1) + dblAbsImg * dblAbsImg); + dblS = dsqrts( (dblAbsReal - 1) * (dblAbsReal - 1) + dblAbsImg * dblAbsImg); + dblA = 0.5 * ( dblR + dblS ); + dblB = dblAbsReal / dblA; + + + /* compute the real part */ + if(dblB <= sdblBcross) + _pdblReal = dasins(dblB); + else if( dblAbsReal <= 1) + _pdblReal = datans(dblAbsReal / dsqrts( 0.5 * (dblA + dblAbsReal) * ( (dblAbsImg * dblAbsImg) / (dblR + (dblAbsReal + 1)) + (dblS + (1 - dblAbsReal))))); + else + _pdblReal = datans(dblAbsReal / (dblAbsImg * dsqrts(0.5 * ((dblA + dblAbsReal) / (dblR + (dblAbsReal + 1)) + (dblA + dblAbsReal) / (dblS + (dblAbsReal-1)))))); + + /* compute the imaginary part */ + if(dblA <= sdblAcross) + { + double dblImg1 = 0; + + if(dblAbsReal < 1) + /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(y**2)/(S+(1.d0-x))) */ + dblImg1 = 0.5 * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + dblAbsImg * dblAbsImg / (dblS + (1 - dblAbsReal))); + else + /* Am1 = 0.5d0*((y**2)/(R+(x+1.d0))+(S+(x-1.d0))) */ + dblImg1 = 0.5 * (dblAbsImg * dblAbsImg / (dblR + (dblAbsReal + 1)) + (dblS + (dblAbsReal - 1))); + /* ai = logp1(Am1 + sqrt(Am1*(A+1.d0))) */ + dblTemp = dblImg1 + dsqrts(dblImg1 * (dblA + 1)); + _pdblImg = dlog1ps(dblTemp); + } + else + /* ai = log(A + sqrt(A**2 - 1.d0)) */ + _pdblImg = dlogs(dblA + dsqrts(dblA * dblA - 1)); + } + else + { + /* evaluation in the special regions ... */ + if(dblAbsImg <= dblEpsm * dabss(dblAbsReal - 1)) + { + if(dblAbsReal < 1) + { + _pdblReal = dasins(dblAbsReal); + _pdblImg = dblAbsImg / dsqrts((1 + dblAbsReal) * (1 - dblAbsReal)); + } + else + { + _pdblReal = sdblPi_2; + if(dblAbsReal <= dblLsup) + { + dblTemp = (dblAbsReal - 1) + dsqrts((dblAbsReal - 1) * (dblAbsReal + 1)); + _pdblImg = dlog1ps(dblTemp); + } + else + _pdblImg = sdblLn2 + dlogs(dblAbsReal); + } + } + else if(dblAbsImg < dblLinf) + { + _pdblReal = sdblPi_2 - dsqrts(dblAbsImg); + _pdblImg = dsqrts(dblAbsImg); + } + else if((dblEpsm * dblAbsImg - 1 >= dblAbsReal)) + { + _pdblReal = dblAbsReal * dblAbsImg; + _pdblImg = sdblLn2 + dlogs(dblAbsReal); + } + else if(dblAbsReal > 1) + { + _pdblReal = datans(dblAbsReal / dblAbsImg); + dblTemp = (dblAbsReal / dblAbsImg) * (dblAbsReal / dblAbsImg); + _pdblImg = sdblLn2 + dlogs(dblAbsReal) + 0.5 * dlog1ps(dblTemp); + } + else + { + double dblTemp2 = dsqrts(1 + dblAbsImg * dblAbsImg); + _pdblReal = dblAbsReal / dblTemp2; + dblTemp = 2 * dblAbsImg * (dblAbsImg + dblTemp2); + _pdblImg = 0.5 * dlog1ps(dblTemp); + } + } + _pdblReal *= iSignReal; + _pdblImg *= iSignImg; + + return (DoubleComplex(_pdblReal, _pdblImg)); +} diff --git a/src/c/elementaryFunctions/atan/catans.c b/src/c/elementaryFunctions/atan/catans.c index d3a232fb..d2081181 100644 --- a/src/c/elementaryFunctions/atan/catans.c +++ b/src/c/elementaryFunctions/atan/catans.c @@ -1,248 +1,249 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -/* - PURPOSE - watan compute the arctangent of a complex number - y = yr + i yi = atan(x), x = xr + i xi - - CALLING LIST / PARAMETERS - subroutine watan(xr,xi,yr,yi) - double precision xr,xi,yr,yi - - xr,xi: real and imaginary parts of the complex number - yr,yi: real and imaginary parts of the result - yr,yi may have the same memory cases than xr et xi - - COPYRIGHT (C) 2001 Bruno Pincon and Lydia van Dijk - Written by Bruno Pincon so - as to get more precision. Also to fix the - behavior at the singular points and at the branch cuts. - Polished by Lydia van Dijk - - - CHANGES : - (Bruno on 2001 May 22) for ysptrk use a - minimax polynome to enlarge the special - evaluation zone |s| < SLIM. Also rename - this function as lnp1m1. - - (Bruno on 2001 June 7) better handling - of spurious over/underflow ; remove - the call to pythag ; better accuracy - in the real part for z near +-i - - EXTERNALS FUNCTIONS - dlamch - lnp1m1 (at the end of this file) - - ALGORITHM : noting z = a + i*b, we have: - Z = yr + yi*b = arctan(z) = (i/2) * log( (i+z)/(i-z) ) - - This function has two branch points at +i and -i and the - chosen branch cuts are the two half-straight lines - D1 = [i, i*oo) and D2 = (-i*oo, i]. The function is then - analytic in C \ (D1 U D2)). - - From the definition it follows that: - - yr = 0.5 Arg ( (i+z)/(i-z) ) (1) - yi = 0.5 log (|(i+z)/(i-z)|) (2) - - so lim (z -> +- i) yr = undefined (and Nan is logical) - lim (z -> +i) yi = +oo - lim (z -> -i) yi = -oo - - The real part of arctan(z) is discontinuous across D1 and D2 - and we impose the following definitions: - if imag(z) > 1 then - Arg(arctan(z)) = pi/2 (=lim real(z) -> 0+) - if imag(z) < 1 then - Arg(arctan(z)) = -pi/2 (=lim real(z) -> 0-) - - - Basic evaluation: if we write (i+z)/(i-z) using - z = a + i*b, we get: - - i+z 1-(a**2+b**2) + i*(2a) - --- = ---------------------- - i-z a**2 + (1-b)**2 - - then, with r2 = |z|^2 = a**2 + b**2 : - - yr = 0.5 * Arg(1-r2 + (2*a)*i) - = 0.5 * atan2(2a, (1-r2)) (3) - - This formula is changed when r2 > RMAX (max pos float) - and also when |1-r2| and |a| are near 0 (see comments - in the code). - - After some math: - - yi = 0.25 * log( (a**2 + (b + 1)**2) / - (a**2 + (b - 1)**2) ) (4) - - Evaluation for "big" |z| - ------------------------ - - If |z| is "big", the direct evaluation of yi by (4) may - suffer of innaccuracies and of spurious overflow. Noting - that s = 2 b / (1 + |z|**2), we have: - - yi = 0.25 log ( (1 + s)/(1 - s) ) (5) - - 3 5 - yi = 0.25*( 2 * ( s + 1/3 s + 1/5 s + ... )) - - yi = 0.25 * lnp1m1(s) if |s| < SLIM - - So if |s| is less than SLIM we switch to a special - evaluation done by the function lnp1m1. The - threshold value SLIM is choosen by experiment - (with the Pari-gp software). For |s| - "very small" we used a truncated taylor dvp, - else a minimax polynome (see lnp1m1). - - To avoid spurious overflows (which result in spurious - underflows for s) in computing s with s= 2 b / (1 + |z|**2) - when |z|^2 > RMAX (max positive float) we use : - - s = 2d0 / ( (a/b)*a + b ) - - but if |b| = Inf this formula leads to NaN when - |a| is also Inf. As we have : - - |s| <= 2 / |b| - - we impose simply : s = 0 when |b| = Inf - - Evaluation for z very near to i or -i: - -------------------------------------- - Floating point numbers of the form a+i or a-i with 0 < - a**2 < tiny (approximately 1d-308) may lead to underflow - (i.e., a**2 = 0) and the logarithm will break formula (4). - So we switch to the following formulas: - - If b = +-1 and |a| < sqrt(tiny) approximately 1d-150 (say) - then (by using that a**2 + 4 = 4 in machine for such a): - - yi = 0.5 * log( 2/|a| ) for b=1 - - yi = 0.5 * log( |a|/2 ) for b=-1 - - finally: yi = 0.5 * sign(b) * log( 2/|a| ) - yi = 0.5 * sign(b) * (log(2) - log(|a|)) (6) - - The last trick is to avoid overflow for |a|=tiny! In fact - this formula may be used until a**2 + 4 = 4 so that the - threshold value may be larger. -*/ - -#include -#include "atan.h" -#include "abs.h" -#include "lnp1m1.h" -#include "lapack.h" - -#define _sign(a, b) b >=0 ? a : -a - -floatComplex catans(floatComplex z) { - static float sSlim = 0.2f; - /* . - ** / \ WARNING : this algorithm was based on double precision - ** / ! \ using float truncate the value to 0. - ** `----' - ** - ** static float sAlim = 1E-150f; - */ - static float sAlim = 0.0f; - static float sTol = 0.3f; - static float sLn2 = 0.6931471805599453094172321f; - - float RMax = (float) getOverflowThreshold(); - float Pi_2 = 2.0f * satans(1); - - float _inReal = creals(z); - float _inImg = cimags(z); - float _outReal = 0; - float _outImg = 0; - - /* Temporary variables */ - float R2 = 0; - float S = 0; - - - if(_inImg == 0) - { - _outReal = satans(_inReal); - _outImg = 0; - } - else - { - R2 = _inReal * _inReal + _inImg * _inImg; /* Oo */ - if(R2 > RMax) - { - if( dabss(_inImg) > RMax) - S = 0; - else - S = 1.0f / (((0.5f * _inReal) / _inImg) * _inReal + 0.5f * _inImg ); - } - else - S = (2 * _inImg) / (1+R2); - - if(dabss(S) < sSlim) - { - /* - s is small: |s| < SLIM <=> |z| outside the following disks: - D+ = D(center = [0; 1/slim], radius = sqrt(1/slim**2 - 1)) if b > 0 - D- = D(center = [0; -1/slim], radius = sqrt(1/slim**2 - 1)) if b < 0 - use the special evaluation of log((1+s)/(1-s)) (5) - */ - _outImg = slnp1m1s(S) * 0.25f; - } - else - { - if(sabss(S) == 1 && sabss(_inReal) <= sAlim) - { - /* |s| >= SLIM => |z| is inside D+ or D- */ - _outImg = _sign(0.5f,_inImg) * ( sLn2 - logf(sabss(_inReal))); - } - else - { - _outImg = 0.25f * logf((powf(_inReal,2) + powf((_inImg + 1.0f),2)) / (powf(_inReal,2) + powf((_inImg - 1.0f),2))); - } - } - if(_inReal == 0) - {/* z is purely imaginary */ - if( dabss(_inImg) > 1) - {/* got sign(b) * pi/2 */ - _outReal = _sign(1, _inImg) * Pi_2; - } - else if( dabss(_inImg) == 1) - {/* got a Nan with 0/0 */ - _outReal = (_inReal - _inReal) / (_inReal - _inReal); /* Oo */ - } - else - _outReal = 0; - } - else if(R2 > RMax) - {/* _outImg is necessarily very near sign(a)* pi/2 */ - _outReal = _sign(1, _inReal) * Pi_2; - } - else if(sabss(1 - R2) + sabss(_inReal) <= sTol) - {/* |b| is very near 1 (and a is near 0) some cancellation occur in the (next) generic formula */ - _outReal = 0.5f * atan2f(2.0f * _inReal, (1.0f - _inImg) * (1.0f + _inImg) - powf(_inReal,2.0f)); - } - else - _outReal = 0.5f * atan2f(2.0f * _inReal, 1.0f - R2); - } - - return FloatComplex(_outReal, _outImg); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +/* + PURPOSE + watan compute the arctangent of a complex number + y = yr + i yi = atan(x), x = xr + i xi + + CALLING LIST / PARAMETERS + subroutine watan(xr,xi,yr,yi) + double precision xr,xi,yr,yi + + xr,xi: real and imaginary parts of the complex number + yr,yi: real and imaginary parts of the result + yr,yi may have the same memory cases than xr et xi + + COPYRIGHT (C) 2001 Bruno Pincon and Lydia van Dijk + Written by Bruno Pincon so + as to get more precision. Also to fix the + behavior at the singular points and at the branch cuts. + Polished by Lydia van Dijk + + + CHANGES : - (Bruno on 2001 May 22) for ysptrk use a + minimax polynome to enlarge the special + evaluation zone |s| < SLIM. Also rename + this function as lnp1m1. + - (Bruno on 2001 June 7) better handling + of spurious over/underflow ; remove + the call to pythag ; better accuracy + in the real part for z near +-i + + EXTERNALS FUNCTIONS + dlamch + lnp1m1 (at the end of this file) + + ALGORITHM : noting z = a + i*b, we have: + Z = yr + yi*b = arctan(z) = (i/2) * log( (i+z)/(i-z) ) + + This function has two branch points at +i and -i and the + chosen branch cuts are the two half-straight lines + D1 = [i, i*oo) and D2 = (-i*oo, i]. The function is then + analytic in C \ (D1 U D2)). + + From the definition it follows that: + + yr = 0.5 Arg ( (i+z)/(i-z) ) (1) + yi = 0.5 log (|(i+z)/(i-z)|) (2) + + so lim (z -> +- i) yr = undefined (and Nan is logical) + lim (z -> +i) yi = +oo + lim (z -> -i) yi = -oo + + The real part of arctan(z) is discontinuous across D1 and D2 + and we impose the following definitions: + if imag(z) > 1 then + Arg(arctan(z)) = pi/2 (=lim real(z) -> 0+) + if imag(z) < 1 then + Arg(arctan(z)) = -pi/2 (=lim real(z) -> 0-) + + + Basic evaluation: if we write (i+z)/(i-z) using + z = a + i*b, we get: + + i+z 1-(a**2+b**2) + i*(2a) + --- = ---------------------- + i-z a**2 + (1-b)**2 + + then, with r2 = |z|^2 = a**2 + b**2 : + + yr = 0.5 * Arg(1-r2 + (2*a)*i) + = 0.5 * atan2(2a, (1-r2)) (3) + + This formula is changed when r2 > RMAX (max pos float) + and also when |1-r2| and |a| are near 0 (see comments + in the code). + + After some math: + + yi = 0.25 * log( (a**2 + (b + 1)**2) / + (a**2 + (b - 1)**2) ) (4) + + Evaluation for "big" |z| + ------------------------ + + If |z| is "big", the direct evaluation of yi by (4) may + suffer of innaccuracies and of spurious overflow. Noting + that s = 2 b / (1 + |z|**2), we have: + + yi = 0.25 log ( (1 + s)/(1 - s) ) (5) + + 3 5 + yi = 0.25*( 2 * ( s + 1/3 s + 1/5 s + ... )) + + yi = 0.25 * lnp1m1(s) if |s| < SLIM + + So if |s| is less than SLIM we switch to a special + evaluation done by the function lnp1m1. The + threshold value SLIM is choosen by experiment + (with the Pari-gp software). For |s| + "very small" we used a truncated taylor dvp, + else a minimax polynome (see lnp1m1). + + To avoid spurious overflows (which result in spurious + underflows for s) in computing s with s= 2 b / (1 + |z|**2) + when |z|^2 > RMAX (max positive float) we use : + + s = 2d0 / ( (a/b)*a + b ) + + but if |b| = Inf this formula leads to NaN when + |a| is also Inf. As we have : + + |s| <= 2 / |b| + + we impose simply : s = 0 when |b| = Inf + + Evaluation for z very near to i or -i: + -------------------------------------- + Floating point numbers of the form a+i or a-i with 0 < + a**2 < tiny (approximately 1d-308) may lead to underflow + (i.e., a**2 = 0) and the logarithm will break formula (4). + So we switch to the following formulas: + + If b = +-1 and |a| < sqrt(tiny) approximately 1d-150 (say) + then (by using that a**2 + 4 = 4 in machine for such a): + + yi = 0.5 * log( 2/|a| ) for b=1 + + yi = 0.5 * log( |a|/2 ) for b=-1 + + finally: yi = 0.5 * sign(b) * log( 2/|a| ) + yi = 0.5 * sign(b) * (log(2) - log(|a|)) (6) + + The last trick is to avoid overflow for |a|=tiny! In fact + this formula may be used until a**2 + 4 = 4 so that the + threshold value may be larger. +*/ + +#include +#include "atan.h" +#include "abs.h" +#include "lnp1m1.h" +#include "lapack.h" + +#define _sign(a, b) b >=0 ? a : -a + +floatComplex catans(floatComplex z) { + static float sSlim = 0.2f; + /* . + ** / \ WARNING : this algorithm was based on double precision + ** / ! \ using float truncate the value to 0. + ** `----' + ** + ** static float sAlim = 1E-150f; + */ + static float sAlim = 0.0f; + static float sTol = 0.3f; + static float sLn2 = 0.6931471805599453094172321f; + + float RMax = (float) getOverflowThreshold(); + float Pi_2 = 2.0f * satans(1); + + float _inReal = creals(z); + float _inImg = cimags(z); + float _outReal = 0; + float _outImg = 0; + + /* Temporary variables */ + float R2 = 0; + float S = 0; + + + if(_inImg == 0) + { + _outReal = satans(_inReal); + _outImg = 0; + } + else + { + R2 = _inReal * _inReal + _inImg * _inImg; /* Oo */ + if(R2 > RMax) + { + if( dabss(_inImg) > RMax) + S = 0; + else + S = 1.0f / (((0.5f * _inReal) / _inImg) * _inReal + 0.5f * _inImg ); + } + else + S = (2 * _inImg) / (1+R2); + + if(dabss(S) < sSlim) + { + /* + s is small: |s| < SLIM <=> |z| outside the following disks: + D+ = D(center = [0; 1/slim], radius = sqrt(1/slim**2 - 1)) if b > 0 + D- = D(center = [0; -1/slim], radius = sqrt(1/slim**2 - 1)) if b < 0 + use the special evaluation of log((1+s)/(1-s)) (5) + */ + _outImg = slnp1m1s(S) * 0.25f; + } + else + { + if(sabss(S) == 1 && sabss(_inReal) <= sAlim) + { + /* |s| >= SLIM => |z| is inside D+ or D- */ + _outImg = _sign(0.5f,_inImg) * ( sLn2 - logf(sabss(_inReal))); + } + else + { + _outImg = 0.25f * logf((powf(_inReal,2) + powf((_inImg + 1.0f),2)) / (powf(_inReal,2) + powf((_inImg - 1.0f),2))); + } + } + if(_inReal == 0) + {/* z is purely imaginary */ + if( dabss(_inImg) > 1) + {/* got sign(b) * pi/2 */ + _outReal = _sign(1, _inImg) * Pi_2; + } + else if( dabss(_inImg) == 1) + {/* got a Nan with 0/0 */ + _outReal = (_inReal - _inReal) / (_inReal - _inReal); /* Oo */ + } + else + _outReal = 0; + } + else if(R2 > RMax) + {/* _outImg is necessarily very near sign(a)* pi/2 */ + _outReal = _sign(1, _inReal) * Pi_2; + } + else if(sabss(1 - R2) + sabss(_inReal) <= sTol) + {/* |b| is very near 1 (and a is near 0) some cancellation occur in the (next) generic formula */ + _outReal = 0.5f * atan2f(2.0f * _inReal, (1.0f - _inImg) * (1.0f + _inImg) - powf(_inReal,2.0f)); + } + else + _outReal = 0.5f * atan2f(2.0f * _inReal, 1.0f - R2); + } + + return FloatComplex(_outReal, _outImg); +} diff --git a/src/c/elementaryFunctions/atan/zatans.c b/src/c/elementaryFunctions/atan/zatans.c index c511d790..4b8e9640 100644 --- a/src/c/elementaryFunctions/atan/zatans.c +++ b/src/c/elementaryFunctions/atan/zatans.c @@ -1,241 +1,242 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -/* - PURPOSE - watan compute the arctangent of a complex number - y = yr + i yi = atan(x), x = xr + i xi - - CALLING LIST / PARAMETERS - subroutine watan(xr,xi,yr,yi) - double precision xr,xi,yr,yi - - xr,xi: real and imaginary parts of the complex number - yr,yi: real and imaginary parts of the result - yr,yi may have the same memory cases than xr et xi - - COPYRIGHT (C) 2001 Bruno Pincon and Lydia van Dijk - Written by Bruno Pincon so - as to get more precision. Also to fix the - behavior at the singular points and at the branch cuts. - Polished by Lydia van Dijk - - - CHANGES : - (Bruno on 2001 May 22) for ysptrk use a - minimax polynome to enlarge the special - evaluation zone |s| < SLIM. Also rename - this function as lnp1m1. - - (Bruno on 2001 June 7) better handling - of spurious over/underflow ; remove - the call to pythag ; better accuracy - in the real part for z near +-i - - EXTERNALS FUNCTIONS - dlamch - lnp1m1 (at the end of this file) - - ALGORITHM : noting z = a + i*b, we have: - Z = yr + yi*b = arctan(z) = (i/2) * log( (i+z)/(i-z) ) - - This function has two branch points at +i and -i and the - chosen branch cuts are the two half-straight lines - D1 = [i, i*oo) and D2 = (-i*oo, i]. The function is then - analytic in C \ (D1 U D2)). - - From the definition it follows that: - - yr = 0.5 Arg ( (i+z)/(i-z) ) (1) - yi = 0.5 log (|(i+z)/(i-z)|) (2) - - so lim (z -> +- i) yr = undefined (and Nan is logical) - lim (z -> +i) yi = +oo - lim (z -> -i) yi = -oo - - The real part of arctan(z) is discontinuous across D1 and D2 - and we impose the following definitions: - if imag(z) > 1 then - Arg(arctan(z)) = pi/2 (=lim real(z) -> 0+) - if imag(z) < 1 then - Arg(arctan(z)) = -pi/2 (=lim real(z) -> 0-) - - - Basic evaluation: if we write (i+z)/(i-z) using - z = a + i*b, we get: - - i+z 1-(a**2+b**2) + i*(2a) - --- = ---------------------- - i-z a**2 + (1-b)**2 - - then, with r2 = |z|^2 = a**2 + b**2 : - - yr = 0.5 * Arg(1-r2 + (2*a)*i) - = 0.5 * atan2(2a, (1-r2)) (3) - - This formula is changed when r2 > RMAX (max pos float) - and also when |1-r2| and |a| are near 0 (see comments - in the code). - - After some math: - - yi = 0.25 * log( (a**2 + (b + 1)**2) / - (a**2 + (b - 1)**2) ) (4) - - Evaluation for "big" |z| - ------------------------ - - If |z| is "big", the direct evaluation of yi by (4) may - suffer of innaccuracies and of spurious overflow. Noting - that s = 2 b / (1 + |z|**2), we have: - - yi = 0.25 log ( (1 + s)/(1 - s) ) (5) - - 3 5 - yi = 0.25*( 2 * ( s + 1/3 s + 1/5 s + ... )) - - yi = 0.25 * lnp1m1(s) if |s| < SLIM - - So if |s| is less than SLIM we switch to a special - evaluation done by the function lnp1m1. The - threshold value SLIM is choosen by experiment - (with the Pari-gp software). For |s| - "very small" we used a truncated taylor dvp, - else a minimax polynome (see lnp1m1). - - To avoid spurious overflows (which result in spurious - underflows for s) in computing s with s= 2 b / (1 + |z|**2) - when |z|^2 > RMAX (max positive float) we use : - - s = 2d0 / ( (a/b)*a + b ) - - but if |b| = Inf this formula leads to NaN when - |a| is also Inf. As we have : - - |s| <= 2 / |b| - - we impose simply : s = 0 when |b| = Inf - - Evaluation for z very near to i or -i: - -------------------------------------- - Floating point numbers of the form a+i or a-i with 0 < - a**2 < tiny (approximately 1d-308) may lead to underflow - (i.e., a**2 = 0) and the logarithm will break formula (4). - So we switch to the following formulas: - - If b = +-1 and |a| < sqrt(tiny) approximately 1d-150 (say) - then (by using that a**2 + 4 = 4 in machine for such a): - - yi = 0.5 * log( 2/|a| ) for b=1 - - yi = 0.5 * log( |a|/2 ) for b=-1 - - finally: yi = 0.5 * sign(b) * log( 2/|a| ) - yi = 0.5 * sign(b) * (log(2) - log(|a|)) (6) - - The last trick is to avoid overflow for |a|=tiny! In fact - this formula may be used until a**2 + 4 = 4 so that the - threshold value may be larger. -*/ - -#include -#include "lapack.h" -#include "atan.h" -#include "abs.h" -#include "lnp1m1.h" - -#define _sign(a, b) b >=0 ? a : -a - -doubleComplex zatans(doubleComplex z) { - static double sSlim = 0.2; - static double sAlim = 1E-150; - static double sTol = 0.3; - static double sLn2 = 0.6931471805599453094172321; - - double RMax = getOverflowThreshold(); - double Pi_2 = 2.0 * datans(1); - - double _inReal = zreals(z); - double _inImg = zimags(z); - double _outReal = 0; - double _outImg = 0; - - /* Temporary variables */ - double R2 = 0; - double S = 0; - - - if(_inImg == 0) - { - _outReal = datans(_inReal); - _outImg = 0; - } - else - { - R2 = _inReal * _inReal + _inImg * _inImg; /* Oo */ - if(R2 > RMax) - { - if( dabss(_inImg) > RMax) - S = 0; - else - S = 1 / (((0.5 * _inReal) / _inImg) * _inReal + 0.5 * _inImg ); - } - else - S = (2 * _inImg) / (1+R2); - - if(dabss(S) < sSlim) - { - /* - s is small: |s| < SLIM <=> |z| outside the following disks: - D+ = D(center = [0; 1/slim], radius = sqrt(1/slim**2 - 1)) if b > 0 - D- = D(center = [0; -1/slim], radius = sqrt(1/slim**2 - 1)) if b < 0 - use the special evaluation of log((1+s)/(1-s)) (5) - */ - _outImg = dlnp1m1s(S) * 0.25; - } - else - { - if(dabss(S) == 1 && dabss(_inReal) <= sAlim) - { - /* |s| >= SLIM => |z| is inside D+ or D- */ - _outImg = _sign(0.5,_inImg) * ( sLn2 - log(dabss(_inReal))); - } - else - { - _outImg = 0.25 * log((pow(_inReal,2) + pow((_inImg + 1),2)) / (pow(_inReal,2) + pow((_inImg - 1),2))); - } - } - if(_inReal == 0) - {/* z is purely imaginary */ - if( dabss(_inImg) > 1) - {/* got sign(b) * pi/2 */ - _outReal = _sign(1, _inImg) * Pi_2; - } - else if( dabss(_inImg) == 1) - {/* got a Nan with 0/0 */ - _outReal = (_inReal - _inReal) / (_inReal - _inReal); /* Oo */ - } - else - _outReal = 0; - } - else if(R2 > RMax) - {/* _outImg is necessarily very near sign(a)* pi/2 */ - _outReal = _sign(1, _inReal) * Pi_2; - } - else if(dabss(1 - R2) + dabss(_inReal) <= sTol) - {/* |b| is very near 1 (and a is near 0) some cancellation occur in the (next) generic formula */ - _outReal = 0.5 * atan2(2 * _inReal, (1-_inImg) * (1 + _inImg) - pow(_inReal,2)); - } - else - _outReal = 0.5 * atan2(2 * _inReal, 1 - R2); - } - - return DoubleComplex(_outReal, _outImg); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +/* + PURPOSE + watan compute the arctangent of a complex number + y = yr + i yi = atan(x), x = xr + i xi + + CALLING LIST / PARAMETERS + subroutine watan(xr,xi,yr,yi) + double precision xr,xi,yr,yi + + xr,xi: real and imaginary parts of the complex number + yr,yi: real and imaginary parts of the result + yr,yi may have the same memory cases than xr et xi + + COPYRIGHT (C) 2001 Bruno Pincon and Lydia van Dijk + Written by Bruno Pincon so + as to get more precision. Also to fix the + behavior at the singular points and at the branch cuts. + Polished by Lydia van Dijk + + + CHANGES : - (Bruno on 2001 May 22) for ysptrk use a + minimax polynome to enlarge the special + evaluation zone |s| < SLIM. Also rename + this function as lnp1m1. + - (Bruno on 2001 June 7) better handling + of spurious over/underflow ; remove + the call to pythag ; better accuracy + in the real part for z near +-i + + EXTERNALS FUNCTIONS + dlamch + lnp1m1 (at the end of this file) + + ALGORITHM : noting z = a + i*b, we have: + Z = yr + yi*b = arctan(z) = (i/2) * log( (i+z)/(i-z) ) + + This function has two branch points at +i and -i and the + chosen branch cuts are the two half-straight lines + D1 = [i, i*oo) and D2 = (-i*oo, i]. The function is then + analytic in C \ (D1 U D2)). + + From the definition it follows that: + + yr = 0.5 Arg ( (i+z)/(i-z) ) (1) + yi = 0.5 log (|(i+z)/(i-z)|) (2) + + so lim (z -> +- i) yr = undefined (and Nan is logical) + lim (z -> +i) yi = +oo + lim (z -> -i) yi = -oo + + The real part of arctan(z) is discontinuous across D1 and D2 + and we impose the following definitions: + if imag(z) > 1 then + Arg(arctan(z)) = pi/2 (=lim real(z) -> 0+) + if imag(z) < 1 then + Arg(arctan(z)) = -pi/2 (=lim real(z) -> 0-) + + + Basic evaluation: if we write (i+z)/(i-z) using + z = a + i*b, we get: + + i+z 1-(a**2+b**2) + i*(2a) + --- = ---------------------- + i-z a**2 + (1-b)**2 + + then, with r2 = |z|^2 = a**2 + b**2 : + + yr = 0.5 * Arg(1-r2 + (2*a)*i) + = 0.5 * atan2(2a, (1-r2)) (3) + + This formula is changed when r2 > RMAX (max pos float) + and also when |1-r2| and |a| are near 0 (see comments + in the code). + + After some math: + + yi = 0.25 * log( (a**2 + (b + 1)**2) / + (a**2 + (b - 1)**2) ) (4) + + Evaluation for "big" |z| + ------------------------ + + If |z| is "big", the direct evaluation of yi by (4) may + suffer of innaccuracies and of spurious overflow. Noting + that s = 2 b / (1 + |z|**2), we have: + + yi = 0.25 log ( (1 + s)/(1 - s) ) (5) + + 3 5 + yi = 0.25*( 2 * ( s + 1/3 s + 1/5 s + ... )) + + yi = 0.25 * lnp1m1(s) if |s| < SLIM + + So if |s| is less than SLIM we switch to a special + evaluation done by the function lnp1m1. The + threshold value SLIM is choosen by experiment + (with the Pari-gp software). For |s| + "very small" we used a truncated taylor dvp, + else a minimax polynome (see lnp1m1). + + To avoid spurious overflows (which result in spurious + underflows for s) in computing s with s= 2 b / (1 + |z|**2) + when |z|^2 > RMAX (max positive float) we use : + + s = 2d0 / ( (a/b)*a + b ) + + but if |b| = Inf this formula leads to NaN when + |a| is also Inf. As we have : + + |s| <= 2 / |b| + + we impose simply : s = 0 when |b| = Inf + + Evaluation for z very near to i or -i: + -------------------------------------- + Floating point numbers of the form a+i or a-i with 0 < + a**2 < tiny (approximately 1d-308) may lead to underflow + (i.e., a**2 = 0) and the logarithm will break formula (4). + So we switch to the following formulas: + + If b = +-1 and |a| < sqrt(tiny) approximately 1d-150 (say) + then (by using that a**2 + 4 = 4 in machine for such a): + + yi = 0.5 * log( 2/|a| ) for b=1 + + yi = 0.5 * log( |a|/2 ) for b=-1 + + finally: yi = 0.5 * sign(b) * log( 2/|a| ) + yi = 0.5 * sign(b) * (log(2) - log(|a|)) (6) + + The last trick is to avoid overflow for |a|=tiny! In fact + this formula may be used until a**2 + 4 = 4 so that the + threshold value may be larger. +*/ + +#include +#include "lapack.h" +#include "atan.h" +#include "abs.h" +#include "lnp1m1.h" + +#define _sign(a, b) b >=0 ? a : -a + +doubleComplex zatans(doubleComplex z) { + static double sSlim = 0.2; + static double sAlim = 1E-150; + static double sTol = 0.3; + static double sLn2 = 0.6931471805599453094172321; + + double RMax = getOverflowThreshold(); + double Pi_2 = 2.0 * datans(1); + + double _inReal = zreals(z); + double _inImg = zimags(z); + double _outReal = 0; + double _outImg = 0; + + /* Temporary variables */ + double R2 = 0; + double S = 0; + + + if(_inImg == 0) + { + _outReal = datans(_inReal); + _outImg = 0; + } + else + { + R2 = _inReal * _inReal + _inImg * _inImg; /* Oo */ + if(R2 > RMax) + { + if( dabss(_inImg) > RMax) + S = 0; + else + S = 1 / (((0.5 * _inReal) / _inImg) * _inReal + 0.5 * _inImg ); + } + else + S = (2 * _inImg) / (1+R2); + + if(dabss(S) < sSlim) + { + /* + s is small: |s| < SLIM <=> |z| outside the following disks: + D+ = D(center = [0; 1/slim], radius = sqrt(1/slim**2 - 1)) if b > 0 + D- = D(center = [0; -1/slim], radius = sqrt(1/slim**2 - 1)) if b < 0 + use the special evaluation of log((1+s)/(1-s)) (5) + */ + _outImg = dlnp1m1s(S) * 0.25; + } + else + { + if(dabss(S) == 1 && dabss(_inReal) <= sAlim) + { + /* |s| >= SLIM => |z| is inside D+ or D- */ + _outImg = _sign(0.5,_inImg) * ( sLn2 - log(dabss(_inReal))); + } + else + { + _outImg = 0.25 * log((pow(_inReal,2) + pow((_inImg + 1),2)) / (pow(_inReal,2) + pow((_inImg - 1),2))); + } + } + if(_inReal == 0) + {/* z is purely imaginary */ + if( dabss(_inImg) > 1) + {/* got sign(b) * pi/2 */ + _outReal = _sign(1, _inImg) * Pi_2; + } + else if( dabss(_inImg) == 1) + {/* got a Nan with 0/0 */ + _outReal = (_inReal - _inReal) / (_inReal - _inReal); /* Oo */ + } + else + _outReal = 0; + } + else if(R2 > RMax) + {/* _outImg is necessarily very near sign(a)* pi/2 */ + _outReal = _sign(1, _inReal) * Pi_2; + } + else if(dabss(1 - R2) + dabss(_inReal) <= sTol) + {/* |b| is very near 1 (and a is near 0) some cancellation occur in the (next) generic formula */ + _outReal = 0.5 * atan2(2 * _inReal, (1-_inImg) * (1 + _inImg) - pow(_inReal,2)); + } + else + _outReal = 0.5 * atan2(2 * _inReal, 1 - R2); + } + + return DoubleComplex(_outReal, _outImg); +} diff --git a/src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c b/src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c index 28cdbe0c..7e1759be 100644 --- a/src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c +++ b/src/c/elementaryFunctions/lnp1m1/dlnp1m1s.c @@ -1,76 +1,77 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -#include "lnp1m1.h" -#include "abs.h" - -/* - PURPOSE : Compute v = log ( (1 + s)/(1 - s) ) - for small s, this is for |s| < SLIM = 0.20 - - ALGORITHM : - 1/ if |s| is "very small" we use a truncated - taylor dvp (by keeping 3 terms) from : - 2 4 6 - t = 2 * s * ( 1 + 1/3 s + 1/5 s + [ 1/7 s + ....] ) - 2 4 - t = 2 * s * ( 1 + 1/3 s + 1/5 s + er) - - The limit E until we use this formula may be simply - gotten so that the negliged part er is such that : - 2 4 - (#) er <= epsm * ( 1 + 1/3 s + 1/5 s ) for all |s|<= E - - As er = 1/7 s^6 + 1/9 s^8 + ... - er <= 1/7 * s^6 ( 1 + s^2 + s^4 + ...) = 1/7 s^6/(1-s^2) - - the inequality (#) is forced if : - - 1/7 s^6 / (1-s^2) <= epsm * ( 1 + 1/3 s^2 + 1/5 s^4 ) - - s^6 <= 7 epsm * (1 - 2/3 s^2 - 3/15 s^4 - 1/5 s^6) - - So that E is very near (7 epsm)^(1/6) (approximately 3.032d-3): - - 2/ For larger |s| we used a minimax polynome : - - yi = s * (2 + d3 s^3 + d5 s^5 .... + d13 s^13 + d15 s^15) - - This polynome was computed (by some remes algorithm) following - (*) the sin(x) example (p 39) of the book : - - "ELEMENTARY FUNCTIONS" - "Algorithms and implementation" - J.M. Muller (Birkhauser) - - (*) without the additionnal raffinement to get the first coefs - very near floating point numbers) -*/ -double dlnp1m1s(double Var) -{ - static double D3 = 0.66666666666672679472; - static double D5 = 0.39999999996176889299; - static double D7 = 0.28571429392829380980; - static double D9 = 0.22222138684562683797; - static double D11 = 0.18186349187499222459; - static double D13 = 0.15250315884469364710; - static double D15 = 0.15367270224757008114; - static double E = 3.032E-3; - static double C3 = 2.0/3.0; - static double C5 = 2.0/5.0; - - double S2 = Var * Var; - if( dabss(Var) <= E) - return Var * (2 + S2 * (C3 + C5 * S2)); - else - return Var * (2 + S2 * (D3 + S2 * (D5 + S2 * (D7 + S2 * (D9 + S2 * (D11 + S2 * (D13 + S2 * D15))))))); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +#include "lnp1m1.h" +#include "abs.h" + +/* + PURPOSE : Compute v = log ( (1 + s)/(1 - s) ) + for small s, this is for |s| < SLIM = 0.20 + + ALGORITHM : + 1/ if |s| is "very small" we use a truncated + taylor dvp (by keeping 3 terms) from : + 2 4 6 + t = 2 * s * ( 1 + 1/3 s + 1/5 s + [ 1/7 s + ....] ) + 2 4 + t = 2 * s * ( 1 + 1/3 s + 1/5 s + er) + + The limit E until we use this formula may be simply + gotten so that the negliged part er is such that : + 2 4 + (#) er <= epsm * ( 1 + 1/3 s + 1/5 s ) for all |s|<= E + + As er = 1/7 s^6 + 1/9 s^8 + ... + er <= 1/7 * s^6 ( 1 + s^2 + s^4 + ...) = 1/7 s^6/(1-s^2) + + the inequality (#) is forced if : + + 1/7 s^6 / (1-s^2) <= epsm * ( 1 + 1/3 s^2 + 1/5 s^4 ) + + s^6 <= 7 epsm * (1 - 2/3 s^2 - 3/15 s^4 - 1/5 s^6) + + So that E is very near (7 epsm)^(1/6) (approximately 3.032d-3): + + 2/ For larger |s| we used a minimax polynome : + + yi = s * (2 + d3 s^3 + d5 s^5 .... + d13 s^13 + d15 s^15) + + This polynome was computed (by some remes algorithm) following + (*) the sin(x) example (p 39) of the book : + + "ELEMENTARY FUNCTIONS" + "Algorithms and implementation" + J.M. Muller (Birkhauser) + + (*) without the additionnal raffinement to get the first coefs + very near floating point numbers) +*/ +double dlnp1m1s(double Var) +{ + static double D3 = 0.66666666666672679472; + static double D5 = 0.39999999996176889299; + static double D7 = 0.28571429392829380980; + static double D9 = 0.22222138684562683797; + static double D11 = 0.18186349187499222459; + static double D13 = 0.15250315884469364710; + static double D15 = 0.15367270224757008114; + static double E = 3.032E-3; + static double C3 = 2.0/3.0; + static double C5 = 2.0/5.0; + + double S2 = Var * Var; + if( dabss(Var) <= E) + return Var * (2 + S2 * (C3 + C5 * S2)); + else + return Var * (2 + S2 * (D3 + S2 * (D5 + S2 * (D7 + S2 * (D9 + S2 * (D11 + S2 * (D13 + S2 * D15))))))); +} diff --git a/src/c/elementaryFunctions/lnp1m1/slnp1m1s.c b/src/c/elementaryFunctions/lnp1m1/slnp1m1s.c index 6c991cc0..9940810c 100644 --- a/src/c/elementaryFunctions/lnp1m1/slnp1m1s.c +++ b/src/c/elementaryFunctions/lnp1m1/slnp1m1s.c @@ -1,76 +1,77 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -#include "lnp1m1.h" -#include "abs.h" - -/* - PURPOSE : Compute v = log ( (1 + s)/(1 - s) ) - for small s, this is for |s| < SLIM = 0.20 - - ALGORITHM : - 1/ if |s| is "very small" we use a truncated - taylor dvp (by keeping 3 terms) from : - 2 4 6 - t = 2 * s * ( 1 + 1/3 s + 1/5 s + [ 1/7 s + ....] ) - 2 4 - t = 2 * s * ( 1 + 1/3 s + 1/5 s + er) - - The limit E until we use this formula may be simply - gotten so that the negliged part er is such that : - 2 4 - (#) er <= epsm * ( 1 + 1/3 s + 1/5 s ) for all |s|<= E - - As er = 1/7 s^6 + 1/9 s^8 + ... - er <= 1/7 * s^6 ( 1 + s^2 + s^4 + ...) = 1/7 s^6/(1-s^2) - - the inequality (#) is forced if : - - 1/7 s^6 / (1-s^2) <= epsm * ( 1 + 1/3 s^2 + 1/5 s^4 ) - - s^6 <= 7 epsm * (1 - 2/3 s^2 - 3/15 s^4 - 1/5 s^6) - - So that E is very near (7 epsm)^(1/6) (approximately 3.032d-3): - - 2/ For larger |s| we used a minimax polynome : - - yi = s * (2 + d3 s^3 + d5 s^5 .... + d13 s^13 + d15 s^15) - - This polynome was computed (by some remes algorithm) following - (*) the sin(x) example (p 39) of the book : - - "ELEMENTARY FUNCTIONS" - "Algorithms and implementation" - J.M. Muller (Birkhauser) - - (*) without the additionnal raffinement to get the first coefs - very near floating point numbers) -*/ -float slnp1m1s(float Var) -{ - static float D3 = 0.66666666666672679472f; - static float D5 = 0.39999999996176889299f; - static float D7 = 0.28571429392829380980f; - static float D9 = 0.22222138684562683797f; - static float D11 = 0.18186349187499222459f; - static float D13 = 0.15250315884469364710f; - static float D15 = 0.15367270224757008114f; - static float E = 3.032E-3f; - static float C3 = 2.0f/3.0f; - static float C5 = 2.0f/5.0f; - - float S2 = Var * Var; - if( sabss(Var) <= E) - return Var * (2 + S2 * (C3 + C5 * S2)); - else - return Var * (2 + S2 * (D3 + S2 * (D5 + S2 * (D7 + S2 * (D9 + S2 * (D11 + S2 * (D13 + S2 * D15))))))); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +#include "lnp1m1.h" +#include "abs.h" + +/* + PURPOSE : Compute v = log ( (1 + s)/(1 - s) ) + for small s, this is for |s| < SLIM = 0.20 + + ALGORITHM : + 1/ if |s| is "very small" we use a truncated + taylor dvp (by keeping 3 terms) from : + 2 4 6 + t = 2 * s * ( 1 + 1/3 s + 1/5 s + [ 1/7 s + ....] ) + 2 4 + t = 2 * s * ( 1 + 1/3 s + 1/5 s + er) + + The limit E until we use this formula may be simply + gotten so that the negliged part er is such that : + 2 4 + (#) er <= epsm * ( 1 + 1/3 s + 1/5 s ) for all |s|<= E + + As er = 1/7 s^6 + 1/9 s^8 + ... + er <= 1/7 * s^6 ( 1 + s^2 + s^4 + ...) = 1/7 s^6/(1-s^2) + + the inequality (#) is forced if : + + 1/7 s^6 / (1-s^2) <= epsm * ( 1 + 1/3 s^2 + 1/5 s^4 ) + + s^6 <= 7 epsm * (1 - 2/3 s^2 - 3/15 s^4 - 1/5 s^6) + + So that E is very near (7 epsm)^(1/6) (approximately 3.032d-3): + + 2/ For larger |s| we used a minimax polynome : + + yi = s * (2 + d3 s^3 + d5 s^5 .... + d13 s^13 + d15 s^15) + + This polynome was computed (by some remes algorithm) following + (*) the sin(x) example (p 39) of the book : + + "ELEMENTARY FUNCTIONS" + "Algorithms and implementation" + J.M. Muller (Birkhauser) + + (*) without the additionnal raffinement to get the first coefs + very near floating point numbers) +*/ +float slnp1m1s(float Var) +{ + static float D3 = 0.66666666666672679472f; + static float D5 = 0.39999999996176889299f; + static float D7 = 0.28571429392829380980f; + static float D9 = 0.22222138684562683797f; + static float D11 = 0.18186349187499222459f; + static float D13 = 0.15250315884469364710f; + static float D15 = 0.15367270224757008114f; + static float E = 3.032E-3f; + static float C3 = 2.0f/3.0f; + static float C5 = 2.0f/5.0f; + + float S2 = Var * Var; + if( sabss(Var) <= E) + return Var * (2 + S2 * (C3 + C5 * S2)); + else + return Var * (2 + S2 * (D3 + S2 * (D5 + S2 * (D7 + S2 * (D9 + S2 * (D11 + S2 * (D13 + S2 * D15))))))); +} diff --git a/src/c/elementaryFunctions/log/clogs.c b/src/c/elementaryFunctions/log/clogs.c index d402d1d1..3b9f6911 100644 --- a/src/c/elementaryFunctions/log/clogs.c +++ b/src/c/elementaryFunctions/log/clogs.c @@ -1,64 +1,65 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -#include -#include "log.h" -#include "lapack.h" -#include "log1p.h" -#include "pythag.h" - -floatComplex clogs(floatComplex in) { - static float sR2 = 1.41421356237309504f; - - float _RealIn = creals(in); - float _ImgIn = cimags(in); - - float _RealOut = 0; - float _ImgOut = 0; - - float RMax = (float) getOverflowThreshold(); - float LInf = sqrtf((float) getUnderflowThreshold()); - float LSup = sqrtf(0.5f * RMax); - - float AbsReal = fabsf(_RealIn); - float AbsImg = fabsf(_ImgIn); - - _ImgOut = atan2f(_ImgIn, _RealIn); - - if(_ImgIn > _RealIn) - {/* switch Real part and Imaginary part */ - float Temp = AbsReal; - AbsReal = AbsImg; - AbsImg = Temp; - } - - if((0.5 <= AbsReal) && (AbsReal <= sR2)) - _RealOut = 0.5f * slog1ps((AbsReal - 1.0f) * (AbsReal + 1.0f) + AbsImg * AbsImg); - else if(LInf < AbsImg && AbsReal < LSup) - _RealOut = 0.5f * slogs(AbsReal * AbsReal + AbsImg * AbsImg); - else if(AbsReal > RMax) - _RealOut = AbsReal; - else - { - float Temp = spythags(AbsReal, AbsImg); - if(Temp <= RMax) - { - _RealOut = slogs(Temp); - } - else /* handle rare spurious overflow with : */ - { - float Temp2 = AbsImg/AbsReal; - _RealOut = slogs(AbsReal) + 0.5f * slog1ps(Temp2 * Temp2); - } - } - return FloatComplex(_RealOut, _ImgOut); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +#include +#include "log.h" +#include "lapack.h" +#include "log1p.h" +#include "pythag.h" + +floatComplex clogs(floatComplex in) { + static float sR2 = 1.41421356237309504f; + + float _RealIn = creals(in); + float _ImgIn = cimags(in); + + float _RealOut = 0; + float _ImgOut = 0; + + float RMax = (float) getOverflowThreshold(); + float LInf = sqrtf((float) getUnderflowThreshold()); + float LSup = sqrtf(0.5f * RMax); + + float AbsReal = fabsf(_RealIn); + float AbsImg = fabsf(_ImgIn); + + _ImgOut = atan2f(_ImgIn, _RealIn); + + if(_ImgIn > _RealIn) + {/* switch Real part and Imaginary part */ + float Temp = AbsReal; + AbsReal = AbsImg; + AbsImg = Temp; + } + + if((0.5 <= AbsReal) && (AbsReal <= sR2)) + _RealOut = 0.5f * slog1ps((AbsReal - 1.0f) * (AbsReal + 1.0f) + AbsImg * AbsImg); + else if(LInf < AbsImg && AbsReal < LSup) + _RealOut = 0.5f * slogs(AbsReal * AbsReal + AbsImg * AbsImg); + else if(AbsReal > RMax) + _RealOut = AbsReal; + else + { + float Temp = spythags(AbsReal, AbsImg); + if(Temp <= RMax) + { + _RealOut = slogs(Temp); + } + else /* handle rare spurious overflow with : */ + { + float Temp2 = AbsImg/AbsReal; + _RealOut = slogs(AbsReal) + 0.5f * slog1ps(Temp2 * Temp2); + } + } + return FloatComplex(_RealOut, _ImgOut); +} diff --git a/src/c/elementaryFunctions/log/zlogs.c b/src/c/elementaryFunctions/log/zlogs.c index dd0792a9..e5e9ded1 100644 --- a/src/c/elementaryFunctions/log/zlogs.c +++ b/src/c/elementaryFunctions/log/zlogs.c @@ -1,64 +1,65 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -#include -#include "log.h" -#include "lapack.h" -#include "log1p.h" -#include "pythag.h" - -doubleComplex zlogs(doubleComplex in) { - static double sR2 = 1.41421356237309504; - - double _RealIn = zreals(in); - double _ImgIn = zimags(in); - - double _RealOut = 0; - double _ImgOut = 0; - - double RMax = getOverflowThreshold(); - double LInf = sqrt(getUnderflowThreshold()); - double LSup = sqrt(0.5 * RMax); - - double AbsReal = fabs(_RealIn); - double AbsImg = fabs(_ImgIn); - - _ImgOut = atan2(_ImgIn, _RealIn); - - if(_ImgIn > _RealIn) - {/* switch Real part and Imaginary part */ - double Temp = AbsReal; - AbsReal = AbsImg; - AbsImg = Temp; - } - - if((0.5 <= AbsReal) && (AbsReal <= sR2)) - _RealOut = 0.5 * dlog1ps((AbsReal - 1) * (AbsReal + 1) + AbsImg * AbsImg); - else if(LInf < AbsImg && AbsReal < LSup) - _RealOut = 0.5 * dlogs(AbsReal * AbsReal + AbsImg * AbsImg); - else if(AbsReal > RMax) - _RealOut = AbsReal; - else - { - double Temp = dpythags(AbsReal, AbsImg); - if(Temp <= RMax) - { - _RealOut = dlogs(Temp); - } - else /* handle rare spurious overflow with : */ - { - double Temp2 = AbsImg/AbsReal; - _RealOut = dlogs(AbsReal) + 0.5 * dlog1ps(Temp2 * Temp2); - } - } - return DoubleComplex(_RealOut, _ImgOut); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2007-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +#include +#include "log.h" +#include "lapack.h" +#include "log1p.h" +#include "pythag.h" + +doubleComplex zlogs(doubleComplex in) { + static double sR2 = 1.41421356237309504; + + double _RealIn = zreals(in); + double _ImgIn = zimags(in); + + double _RealOut = 0; + double _ImgOut = 0; + + double RMax = getOverflowThreshold(); + double LInf = sqrt(getUnderflowThreshold()); + double LSup = sqrt(0.5 * RMax); + + double AbsReal = fabs(_RealIn); + double AbsImg = fabs(_ImgIn); + + _ImgOut = atan2(_ImgIn, _RealIn); + + if(_ImgIn > _RealIn) + {/* switch Real part and Imaginary part */ + double Temp = AbsReal; + AbsReal = AbsImg; + AbsImg = Temp; + } + + if((0.5 <= AbsReal) && (AbsReal <= sR2)) + _RealOut = 0.5 * dlog1ps((AbsReal - 1) * (AbsReal + 1) + AbsImg * AbsImg); + else if(LInf < AbsImg && AbsReal < LSup) + _RealOut = 0.5 * dlogs(AbsReal * AbsReal + AbsImg * AbsImg); + else if(AbsReal > RMax) + _RealOut = AbsReal; + else + { + double Temp = dpythags(AbsReal, AbsImg); + if(Temp <= RMax) + { + _RealOut = dlogs(Temp); + } + else /* handle rare spurious overflow with : */ + { + double Temp2 = AbsImg/AbsReal; + _RealOut = dlogs(AbsReal) + 0.5 * dlog1ps(Temp2 * Temp2); + } + } + return DoubleComplex(_RealOut, _ImgOut); +} diff --git a/src/c/elementaryFunctions/log1p/dlog1ps.c b/src/c/elementaryFunctions/log1p/dlog1ps.c index 4c055c86..e75a05fd 100644 --- a/src/c/elementaryFunctions/log1p/dlog1ps.c +++ b/src/c/elementaryFunctions/log1p/dlog1ps.c @@ -1,33 +1,34 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -#include "log1p.h" -#include "log.h" -#include "lnp1m1.h" - -double dlog1ps(double in) { - static double A = -1.0/3.0; - static double B = 0.5; - - if(in < -1) - {/* got NaN */ - return (in - in) / (in - in); /* NaN */ - } - else if(A <= in && in <= B) - {/* use the function log((1+g)/(1-g)) with g = x/(x + 2) */ - return dlnp1m1s(in / (in + 2)); - } - else - {/* use the standard formula */ - return dlogs(in + 1); - } -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +#include "log1p.h" +#include "log.h" +#include "lnp1m1.h" + +double dlog1ps(double in) { + static double A = -1.0/3.0; + static double B = 0.5; + + if(in < -1) + {/* got NaN */ + return (in - in) / (in - in); /* NaN */ + } + else if(A <= in && in <= B) + {/* use the function log((1+g)/(1-g)) with g = x/(x + 2) */ + return dlnp1m1s(in / (in + 2)); + } + else + {/* use the standard formula */ + return dlogs(in + 1); + } +} diff --git a/src/c/elementaryFunctions/log1p/slog1ps.c b/src/c/elementaryFunctions/log1p/slog1ps.c index 501de6b0..04786524 100644 --- a/src/c/elementaryFunctions/log1p/slog1ps.c +++ b/src/c/elementaryFunctions/log1p/slog1ps.c @@ -1,33 +1,35 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -#include "log1p.h" -#include "log.h" -#include "lnp1m1.h" - -float slog1ps(float in) { - static double A = -1.0/3.0; - static double B = 0.5; - - if(in < -1) - {/* got NaN */ - return (in - in) / (in - in); /* NaN */ - } - else if(A <= in && in <= B) - {/* use the function log((1+g)/(1-g)) with g = x/(x + 2) */ - return slnp1m1s(in / (in + 2)); - } - else - {/* use the standard formula */ - return slogs(in + 1); - } -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +#include "log1p.h" +#include "log.h" +#include "lnp1m1.h" + +float slog1ps(float in) { + static double A = -1.0/3.0; + static double B = 0.5; + + if(in < -1) + {/* got NaN */ + return (in - in) / (in - in); /* NaN */ + } + else if(A <= in && in <= B) + {/* use the function log((1+g)/(1-g)) with g = x/(x + 2) */ + return slnp1m1s(in / (in + 2)); + } + else + {/* use the standard formula */ + return slogs(in + 1); + } +} diff --git a/src/c/elementaryFunctions/sqrt/csqrts.c b/src/c/elementaryFunctions/sqrt/csqrts.c index a2ed819e..a24f9558 100644 --- a/src/c/elementaryFunctions/sqrt/csqrts.c +++ b/src/c/elementaryFunctions/sqrt/csqrts.c @@ -1,6 +1,7 @@ /* * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon * * This file must be used under the terms of the CeCILL. * This source file is licensed as described in the file COPYING, which diff --git a/src/c/elementaryFunctions/sqrt/zsqrts.c b/src/c/elementaryFunctions/sqrt/zsqrts.c index f1a2b05c..3637ddd6 100644 --- a/src/c/elementaryFunctions/sqrt/zsqrts.c +++ b/src/c/elementaryFunctions/sqrt/zsqrts.c @@ -1,6 +1,7 @@ /* * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab * Copyright (C) 2008-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon * * This file must be used under the terms of the CeCILL. * This source file is licensed as described in the file COPYING, which diff --git a/src/c/elementaryFunctions/tan/ztans.c b/src/c/elementaryFunctions/tan/ztans.c index ff130dcc..761da36b 100644 --- a/src/c/elementaryFunctions/tan/ztans.c +++ b/src/c/elementaryFunctions/tan/ztans.c @@ -1,103 +1,104 @@ -/* - * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab - * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET - * - * This file must be used under the terms of the CeCILL. - * This source file is licensed as described in the file COPYING, which - * you should have received as part of this distribution. The terms - * are also available at - * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt - * - */ - -/* - ALGORITHM - based on the formula : - - 0.5 sin(2 xr) + i 0.5 sinh(2 xi) - tan(xr + i xi) = --------------------------------- - cos(xr)^2 + sinh(xi)^2 - - noting d = cos(xr)^2 + sinh(xi)^2, we have : - - yr = 0.5 * sin(2 * xr) / d (1) - - yi = 0.5 * sinh(2 * xi) / d (2) - - to avoid spurious overflows in computing yi with - formula (2) (which results in NaN for yi) - we use also the following formula : - - yi = sign(xi) when |xi| > LIM (3) - - Explanations for (3) : - - we have d = sinh(xi)^2 ( 1 + (cos(xr)/sinh(xi))^2 ), - so when : - - (cos(xr)/sinh(xi))^2 < epsm ( epsm = max relative error - for coding a real in a f.p. - number set F(b,p,emin,emax) - epsm = 0.5 b^(1-p) ) - which is forced when : - - 1/sinh(xi)^2 < epsm (4) - <=> |xi| > asinh(1/sqrt(epsm)) (= 19.06... in ieee 754 double) - - sinh(xi)^2 is a good approximation for d (relative to the f.p. - arithmetic used) and then yr may be approximate with : - - yr = cosh(xi)/sinh(xi) - = sign(xi) (1 + exp(-2 |xi|))/(1 - exp(-2|xi|)) - = sign(xi) (1 + 2 u + 2 u^2 + 2 u^3 + ...) - - with u = exp(-2 |xi|)). Now when : - - 2 exp(-2|xi|) < epsm (2) - <=> |xi| > 0.5 * log(2/epsm) (= 18.71... in ieee 754 double) - - sign(xi) is a good approximation for yr. - - Constraint (1) is stronger than (2) and we take finaly - - LIM = 1 + log(2/sqrt(epsm)) - - (log(2/sqrt(epsm)) being very near asinh(1/sqrt(epsm)) - -AUTHOR - Bruno Pincon -*/ - -#include -#include "lapack.h" -#include "tan.h" -#include "sqrt.h" -#include "cos.h" -#include "sinh.h" -#include "sin.h" -#include "log.h" -#include "abs.h" - -#define localSign(x) x >= 0 ? 1.0 : -1.0 - -doubleComplex ztans(doubleComplex z) { - double Temp = 0; - double Lim = 1 + dlogs(2.0 / dsqrts( getRelativeMachinePrecision())); - double RealIn = zreals(z); - double ImagIn = zimags(z); - double RealOut = 0; - double ImagOut = 0; - - Temp = pow(dcoss(RealIn), 2) + pow(dsinhs(ImagIn), 2); - RealOut = 0.5 * dsins(2 * RealIn) / Temp; - if(dabss(ImagIn) < Lim) - { - ImagOut = 0.5 * dsinhs(2 * ImagIn) / Temp; - } - else - { - ImagOut = localSign(ImagIn); - } - - return DoubleComplex(RealOut, ImagOut); -} +/* + * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab + * Copyright (C) 2006-2008 - INRIA - Bruno JOFRET + * Copyright (C) Bruno Pincon + * + * This file must be used under the terms of the CeCILL. + * This source file is licensed as described in the file COPYING, which + * you should have received as part of this distribution. The terms + * are also available at + * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + * + */ + +/* + ALGORITHM + based on the formula : + + 0.5 sin(2 xr) + i 0.5 sinh(2 xi) + tan(xr + i xi) = --------------------------------- + cos(xr)^2 + sinh(xi)^2 + + noting d = cos(xr)^2 + sinh(xi)^2, we have : + + yr = 0.5 * sin(2 * xr) / d (1) + + yi = 0.5 * sinh(2 * xi) / d (2) + + to avoid spurious overflows in computing yi with + formula (2) (which results in NaN for yi) + we use also the following formula : + + yi = sign(xi) when |xi| > LIM (3) + + Explanations for (3) : + + we have d = sinh(xi)^2 ( 1 + (cos(xr)/sinh(xi))^2 ), + so when : + + (cos(xr)/sinh(xi))^2 < epsm ( epsm = max relative error + for coding a real in a f.p. + number set F(b,p,emin,emax) + epsm = 0.5 b^(1-p) ) + which is forced when : + + 1/sinh(xi)^2 < epsm (4) + <=> |xi| > asinh(1/sqrt(epsm)) (= 19.06... in ieee 754 double) + + sinh(xi)^2 is a good approximation for d (relative to the f.p. + arithmetic used) and then yr may be approximate with : + + yr = cosh(xi)/sinh(xi) + = sign(xi) (1 + exp(-2 |xi|))/(1 - exp(-2|xi|)) + = sign(xi) (1 + 2 u + 2 u^2 + 2 u^3 + ...) + + with u = exp(-2 |xi|)). Now when : + + 2 exp(-2|xi|) < epsm (2) + <=> |xi| > 0.5 * log(2/epsm) (= 18.71... in ieee 754 double) + + sign(xi) is a good approximation for yr. + + Constraint (1) is stronger than (2) and we take finaly + + LIM = 1 + log(2/sqrt(epsm)) + + (log(2/sqrt(epsm)) being very near asinh(1/sqrt(epsm)) + +AUTHOR + Bruno Pincon +*/ + +#include +#include "lapack.h" +#include "tan.h" +#include "sqrt.h" +#include "cos.h" +#include "sinh.h" +#include "sin.h" +#include "log.h" +#include "abs.h" + +#define localSign(x) x >= 0 ? 1.0 : -1.0 + +doubleComplex ztans(doubleComplex z) { + double Temp = 0; + double Lim = 1 + dlogs(2.0 / dsqrts( getRelativeMachinePrecision())); + double RealIn = zreals(z); + double ImagIn = zimags(z); + double RealOut = 0; + double ImagOut = 0; + + Temp = pow(dcoss(RealIn), 2) + pow(dsinhs(ImagIn), 2); + RealOut = 0.5 * dsins(2 * RealIn) / Temp; + if(dabss(ImagIn) < Lim) + { + ImagOut = 0.5 * dsinhs(2 * ImagIn) / Temp; + } + else + { + ImagOut = localSign(ImagIn); + } + + return DoubleComplex(RealOut, ImagOut); +} -- cgit