Cider Isn't Darwin Emulation, Really
Something went wrong. Try again.
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853/* * Copyright (c) 2002 Apple Computer, Inc. All rights reserved. * * @APPLE_LICENSE_HEADER_START@ * * The contents of this file constitute Original Code as defined in and * are subject to the Apple Public Source License Version 1.1 (the * "License"). You may not use this file except in compliance with the * License. Please obtain a copy of the License at * http://www.apple.com/publicsource and read it before using this file. * * This Original Code and all software distributed under the License are * distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, EITHER * EXPRESS OR IMPLIED, AND APPLE HEREBY DISCLAIMS ALL SUCH WARRANTIES, * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE OR NON-INFRINGEMENT. Please see the * License for the specific language governing rights and limitations * under the License. * * @APPLE_LICENSE_HEADER_END@ *//****************************************************************************** File: complex.c**** Contains: C source code for implementations of floating-point** (double) complex functions defined in header file** "complex.h" for PowerPC Macintoshes in native mode.** Transcendental function algorithms are based on the** paper "Branch Cuts for Complex Elementary Functions"** by W. Kahan, May 17, 1987, and on Pascal and C source** code for the SANE 80-/96-bit extended type by Kenton** Hanson and Paul Finlayson, respectively.**** ** Written by: Jon Okada, SANEitation Engineer, ext. 4-4838** ** Copyright: c 1987-1993 by Apple Computer, Inc., all rights reserved.** ** Change History (most recent first):**** 25 Aug 93 ali Changed clog to cLog to avoid clashing with the** stream i/o definition clog.** 14 Jul 93 ali Added #pragma fenv_access on** 22 Feb 93 ali Added a nomaf #pragma.** 05 Feb 93 JPO Modified calls to feclearexcept and feraiseexcept** to reflect changes in "fenv.h".** 18 Dec 92 JPO First created.** ****************************************************************************/
#include "math.h"#include "complex.h"#include "fenv.h"#include "xmmLibm_prefix.h"
#define Real(z) (__real__ z)#define Imag(z) (__imag__ z)
/**************************************************************************** CONSTANTS used by complex functions
#include <stdio.h>#include <math.h>#include <float.h>main(){
float FPKASINHOM4f = asinhf(nextafterf(INFINITY,0.0f))/4.0f;float FPKTHETAf = sqrtf(nextafterf(INFINITY,0.0f))/4.0f;float FPKRHOf = 1.0f/FPKTHETAf;float FPKLOVEREf = FLT_MIN/FLT_EPSILON;
printf("FPKASINHOM4 %16.7e %x\n", FPKASINHOM4f, *(int *)(&FPKASINHOM4f));printf("FPKTHETA %16.7e %x\n", FPKTHETAf, *(int *)(&FPKTHETAf));printf("FPKRHO %16.7e %x\n", FPKRHOf, *(int *)(&FPKRHOf));printf("FPKLOVERE %16.7e %x\n", FPKLOVEREf, *(int *)(&FPKLOVEREf));} static const // underflow threshold / round threshold hexdouble FPKLOVERE = HEXDOUBLE(0x03600000,0x00000000);
static const // underflow threshold / round threshold hexsingle FPKLOVEREf = { 0xc000000 };
static const // exp(709.0) hexdouble FPKEXP709 = HEXDOUBLE(0x7fdd422d,0x2be5dc9b);
static const // asinh(nextafterd(+infinity,0.0))/4.0 hexdouble FPKASINHOM4 = HEXDOUBLE(0x406633ce,0x8fb9f87e);
static const // asinh(nextafterf(+infinity,0.0))/4.0 hexsingle FPKASINHOM4f = { 0x41b2d4fc };
static const // sqrt(nextafterd(+infinity,0.0))/4.0 hexdouble FPKTHETA = HEXDOUBLE(0x5fcfffff,0xffffffff);
static const // sqrt(nextafterf(+infinity,0.0))/4.0 hexsingle FPKTHETAf = { 0x5e7fffff };
static const // 4.0/sqrt(nextafterd(+infinity,0.0)) hexdouble FPKRHO = HEXDOUBLE(0x20100000,0x00000000);
static const // 4.0/sqrt(nextafterf(+infinity,0.0)) hexsingle FPKRHOf = { 0x20800001 };
****************************************************************************/
static const double expOverflowThreshold_d = 0x1.62e42fefa39efp+9;static const double expOverflowValue_d = 0x1.fffffffffff2ap+1023; // exp(overflowThreshold)static const double twiceExpOverflowThresh_d = 0x1.62e42fefa39efp+10;
static const long double expOverflowThreshold_ld = 0xb.17217f7d1cf79abp+10L;static const long double expOverflowValue_ld = 0xf.fffffffffffcd87p+16380L; // exp(overflowThreshold)static const long double twiceExpOverflowThresh_ld = 0xb.17217f7d1cf79abp+10L;
static const double FPKASINHOM4 = 0x1.633ce8fb9f87ep+7;static const float FPKASINHOM4f = 0x1.65a9f8p+4f;static const double FPKTHETA = 0x1.fffffffffffffp+509;static const float FPKTHETAf = 0x1.fffffep+61f;static const double FPKRHO = 0x1p-510;static const float FPKRHOf = 0x1.000002p-62f;
staticdouble complex xdivc( double x, double complex y ) /* returns (real x) / (complex y) */{ double complex z; double r, denom; if ( __builtin_fabs(Real(y)) >= __builtin_fabs(Imag(y)) ) { /* |Real(y)| >= |Imag(y)| */ if (__builtin_fabs(Real(y)) == INFINITY) { /* Imag(y) and Real(y) are infinite */ Real(z) = __builtin_copysign(0.0,Real(y)); Imag(z) = __builtin_copysign(0.0,-Imag(y)); } else { /* |Real(y)| >= finite |Imag(y)| */ r = Imag(y)/Real(y); denom = Real(y) + Imag(y)*r; Real(z) = x/denom; Imag(z) = (-x*r)/denom; } } else { /* |Real(y)| !>= |Imag(y)| */ r = Real(y)/Imag(y); denom = r*Real(y) + Imag(y); Real(z) = (r*x)/denom; Imag(z) = -x/denom; } return z;}
staticfloat complex xdivcf( float x, float complex y ) /* returns (real x) / (complex y) */{ float complex z; float r, denom; if ( __builtin_fabsf(Real(y)) >= __builtin_fabsf(Imag(y)) ) { /* |Real(y)| >= |Imag(y)| */ if (__builtin_fabsf(Real(y)) == INFINITY) { /* Imag(y) and Real(y) are infinite */ Real(z) = __builtin_copysignf(0.0f,Real(y)); Imag(z) = __builtin_copysignf(0.0f,-Imag(y)); } else { /* |Real(y)| >= finite |Imag(y)| */ r = Imag(y)/Real(y); denom = Real(y) + Imag(y)*r; Real(z) = x/denom; Imag(z) = (-x*r)/denom; } } else { /* |Real(y)| !>= |Imag(y)| */ r = Real(y)/Imag(y); denom = r*Real(y) + Imag(y); Real(z) = (r*x)/denom; Imag(z) = -x/denom; } return z;}
staticlong double complex xdivcl( long double x, long double complex y ) /* returns (real x) / (complex y) */{ long double complex z; long double r, denom; if ( __builtin_fabsl(Real(y)) >= __builtin_fabsl(Imag(y)) ) { /* |Real(y)| >= |Imag(y)| */ if (__builtin_fabsl(Real(y)) == INFINITY) { /* Imag(y) and Real(y) are infinite */ Real(z) = __builtin_copysignl(0.0L,Real(y)); Imag(z) = __builtin_copysignl(0.0L,-Imag(y)); } else { /* |Real(y)| >= finite |Imag(y)| */ r = Imag(y)/Real(y); denom = Real(y) + Imag(y)*r; Real(z) = x/denom; Imag(z) = (-x*r)/denom; } } else { /* |Real(y)| !>= |Imag(y)| */ r = Real(y)/Imag(y); denom = r*Real(y) + Imag(y); Real(z) = (r*x)/denom; Imag(z) = -x/denom; } return z;}
/**************************************************************************** double cabs(double complex z) returns the absolute value (magnitude) of its complex argument z, avoiding spurious overflow, underflow, and invalid exceptions. The code is identical to hypot[fl]. On Intel, the cabs functions reside in hypot[fl].s****************************************************************************/
// PowerPC implementation of cabs is here in the ppc complex.c file
/**************************************************************************** double carg(double complex z) returns the argument (in radians) of its complex argument z. The algorithm is from Kahan's paper. The argument of a complex number z = x + i*y is defined as atan2(y,x) for finite x and y. CONSTANTS: FPKPI2 = pi/2.0 to double precision FPKPI = pi to double precision Calls: fpclassify, copysign, fabs, atan****************************************************************************/
double carg ( double complex z ) { return atan2(Imag(z), Real(z)); } float cargf ( float complex z ) { return atan2f(Imag(z), Real(z)); }
long double cargl ( long double complex z ) { return atan2l(Imag(z), Real(z)); }
/**************************************************************************** double complex csqrt(double complex z) returns the complex square root of its argument. The algorithm, which is from the Kahan paper, uses the following identities: sqrt(x + i*y) = sqrt((|z| + Real(z))/2) + i*sqrt((|z| - Real(z))/2) and sqrt(x - i*y) = sqrt((|z| + Real(z))/2) - i*sqrt((|z| - Real(z))/2), where y is positive and x may be positive or negative. CONSTANTS: FPKINF = +infinity Calls: cssqs, scalb, fabs, sqrt, copysign.****************************************************************************/
/* New Intel code written 9/26/06 by scanon * Uses extra precision to compute |z| instead of cssqs(), saving environment calls. * Note that we could also rescale using bits of the SSE2 code from ian's original intel hypot() routine, and will probably * want to do exactly that in the future, to move away from using x87 for this. */
double complex csqrt ( double complex z ){ static const double inf = __builtin_inf(); double u,v; // Special case for infinite y: if (__builtin_fabs(Imag(z)) == inf) return inf + I*Imag(z); // csqrt(x � i�) = � � i� for all x, including NaN. // Special cases for y = NaN: if (Imag(z) != Imag(z)) { if (Real(z) != Real(z)) // csqrt(NaN + iNaN) = NaN + iNaN, quietly. return z; else if (Real(z) == inf) // csqrt(� + iNaN) = � + iNaN return z; else if (Real(z) == -inf) // csqrt(-� � iNaN) = NaN � i�. return Imag(z) + I*__builtin_copysign(inf,Imag(z)); else { // csqrt(x + iNaN) = NaN + iNaN if x is finite. return Imag(z) + I*Imag(z); } } // At this point, we know that y is finite. Deal with special cases for x: // Special case for x = NaN: if (Real(z) != Real(z)) { // csqrt(NaN + ix) = NaN + iNaN. return Real(z) + I*__builtin_copysign(Real(z),Imag(z)); } // Special cases for x = 0: if (Real(z) == 0.0) { if (Imag(z) == 0.0) // csqrt(�0 + i0) = 0 + i0, csqrt(�0 - i0) = 0 - i0. return I*Imag(z); else { // csqrt(0 � iy) = sqrt(y/2) � i sqrt(y/2). u = __builtin_sqrt(0.5*__builtin_fabs(Imag(z))); return u + I*__builtin_copysign(u, Imag(z) ); } } // Special cases for infinte x: if (Real(z) == inf) // csqrt(� � iy) = � � i0 for finite y. return inf + I*__builtin_copysign(0.0,Imag(z)); if (Real(z) == -inf) // csqrt(-� � iy) = 0 � i� for finite y. return I*__builtin_copysign(inf,Imag(z)); // At this point, we know that x is finite, non-zero and y is finite. else { // We use extended (80-bit) precision to avoid over- or under-flow in computing |z|. long double x = __builtin_fabsl(Real(z)); long double y = Imag(z); /* Compute * +---------------- +---------------- * | |z| + |Re z| Im z | |z| - |Re z| * u = | -------------- v = ------ = � | -------------- * \| 2 2u \| 2 * * There is no risk of drastic loss of precision due to cancellation using these formulas, * as there would be if we used the second expression (involving the square root) for v. * * Overflow or Underflow is possible, but only if the actual result does not fit in double width. */ u = (double)__builtin_sqrtl(0.5L*(__builtin_sqrtl(x*x + y*y) + x)); v = 0.5 * (Imag(z) / u); /* If x < 0, then sqrt(z) = |v| + I*copysign(u, Im z). * Otherwise, sqrt(z) = u + I*v. */ if (Real(z) < 0.0) { return __builtin_fabs(v) + I*__builtin_copysign(u,Imag(z)); } else { return u + I*v; } }}
float complex csqrtf ( float complex z ) { static const float inf = __builtin_inff(); float u,v; // Special case for infinite y: if (__builtin_fabsf(Imag(z)) == inf) return inf + I*Imag(z); // csqrt(x � i�) = � � i� for all x, including NaN. // Special cases for y = NaN: if (Imag(z) != Imag(z)) { if (Real(z) != Real(z)) // csqrt(NaN + iNaN) = NaN + iNaN, quietly. return z; else if (Real(z) == inf) // csqrt(� + iNaN) = � + iNaN return z; else if (Real(z) == -inf) // csqrt(-� � iNaN) = NaN � i�. return Imag(z) + I*__builtin_copysignf(inf,Imag(z)); else { // csqrt(x + iNaN) = NaN + iNaN if x is finite. return Imag(z) + I*Imag(z); } } // At this point, we know that y is finite. Deal with special cases for x: // Special case for x = NaN: if (Real(z) != Real(z)) { // csqrt(NaN + ix) = NaN + iNaN. return Real(z) + I*__builtin_copysignf(Real(z),Imag(z)); } // Special cases for x = 0: if (Real(z) == 0.0f) { if (Imag(z) == 0.0f) // csqrt(�0 + i0) = 0 + i0, csqrt(�0 - i0) = 0 - i0. return I*Imag(z); else { // csqrt(0 � iy) = sqrt(y/2) � i sqrt(y/2). u = __builtin_sqrtf(0.5f*__builtin_fabsf(Imag(z))); return u + I*__builtin_copysignf(u, Imag(z) ); } } // Special cases for infinte x: if (Real(z) == inf) // csqrt(� � iy) = � � i0 for finite y. return inf + I*__builtin_copysignf(0.0f,Imag(z)); if (Real(z) == -inf) // csqrt(-� � iy) = 0 � i� for finite y. return I*__builtin_copysignf(inf,Imag(z)); // At this point, we know that x is finite, non-zero and y is finite. else { // We use double (64-bit) precision to avoid over- or under-flow in computing |z|. double x = __builtin_fabs(Real(z)); double y = Imag(z); /* Compute * +---------------- +---------------- * | |z| + |Re z| Im z | |z| - |Re z| * u = | -------------- v = ------ = � | -------------- * \| 2 2u \| 2 * * There is no risk of drastic loss of precision due to cancellation using these formulas, * as there would be if we used the second expression (involving the square root) for v. * * Overflow or Underflow is possible, but only if the actual result does not fit in double width. */ u = (float)__builtin_sqrt(0.5*(__builtin_sqrt(x*x + y*y) + x)); v = 0.5f * (Imag(z) / u); /* If x < 0, then sqrt(z) = |v| + I*copysign(u, Im z). * Otherwise, sqrt(z) = u + I*v. */ if (Real(z) < 0.0f) { return __builtin_fabsf(v) + I*__builtin_copysignf(u,Imag(z)); } else { return u + I*v; } } }
typedef union{ long double ld; struct { uint64_t mantissa; int16_t sexp; };}ld_parts;
long double complex csqrtl ( long double complex z ) { static const long double inf = __builtin_infl(); static const long double zero = 0.0l; static const long double half = 0.5l; long double u,v; // Special case for infinite y: if (__builtin_fabsl(Imag(z)) == inf) return inf + I*Imag(z); // csqrt(x � i�) = � � i� for all x, including NaN. // Special cases for y = NaN: if (Imag(z) != Imag(z)) { if (Real(z) != Real(z)) // csqrt(NaN + iNaN) = NaN + iNaN, quietly. return z; else if (Real(z) == inf) // csqrt(� + iNaN) = � + iNaN return z; else if (Real(z) == -inf) // csqrt(-� � iNaN) = NaN � i�. return Imag(z) + I*__builtin_copysignl(inf,Imag(z)); else { // csqrt(x + iNaN) = NaN + iNaN if x is finite. return Imag(z) + I*Imag(z); } } // Special cases for y = 0: if (Imag(z) == zero) { if (Real(z) == zero) // csqrt(�0 + i0) = 0 + i0, csqrt(�0 - i0) = 0 - i0. return I*Imag(z); else { u = __builtin_sqrtl(__builtin_fabsl(Real(z))); if (Real(z) < zero) return zero + I*__builtin_copysignl(u,Imag(z)); else return u + I*__builtin_copysignl(zero,Imag(z)); } } // At this point, we know that y is finite. Deal with special cases for x: // Special case for x = NaN: if (Real(z) != Real(z)) { // csqrt(NaN + ix) = NaN + iNaN. return Real(z) + I*__builtin_copysignl(Real(z),Imag(z)); } // Special cases for x = 0: if (Real(z) == zero) { // csqrt(0 � iy) = sqrt(y/2) � i sqrt(y/2). u = __builtin_sqrtl(half*__builtin_fabsl(Imag(z))); return u + I*__builtin_copysignl(u, Imag(z) ); } // Special cases for infinte x: if (Real(z) == inf) // csqrt(� � iy) = � � i0 for finite y. return inf + I*__builtin_copysignl(zero,Imag(z)); if (Real(z) == -inf) // csqrt(-� � iy) = 0 � i� for finite y. return I*__builtin_copysignl(inf,Imag(z)); // At this point, we know that x and y are finite, non-zero. else { long double x = __builtin_fabsl(Real(z)); long double y = __builtin_fabsl(Imag(z)); /* Compute * +---------------- +---------------- * | |z| + |Re z| Im z | |z| - |Re z| * u = | -------------- v = ------ = � | -------------- * \| 2 2u \| 2 * * There is no risk of drastic loss of precision due to cancellation using these formulas, * as there would be if we used the second expression (involving the square root) for v. * */ // Scaling code taken from hypotl pretty much wholesale. ld_parts *large = (ld_parts*) &x; ld_parts *small = (ld_parts*) &y; if (large->ld < small->ld) { ld_parts *p = large; large = small; small = p; } int lexp = large->sexp; int sexp = small->sexp; if( lexp == 0 ) { large->ld = large->mantissa; lexp = large->sexp - 16445; } if( sexp == 0 ) { small->ld = small->mantissa; sexp = small->sexp - 16445; } large->sexp = 0x3fff; int scale = 0x3fff - lexp; int small_scale = sexp + scale; if( small_scale < 64 ) small_scale = 64; small->sexp = small_scale; u = __builtin_sqrtl( large->ld * large->ld + small->ld * small->ld ) + x; if (scale%2) scale = 0x3fff - (scale + 1)/2; else { scale = 0x3fff - (scale/2 + 1); u = u + u; } u = __builtin_sqrtl(u); // Rescale result. large->sexp = scale; large->mantissa = 0x8000000000000000ULL; u *= large->ld; // End scaling code. At this point u = sqrt((|z| + |Re z|) / 2). v = Imag(z) / (2.0l * u); /* If x < 0, then sqrt(z) = |v| + I*copysign(u, Im z). * Otherwise, sqrt(z) = u + I*v. */ if (Real(z) < zero) { return __builtin_fabsl(v) + I*__builtin_copysignl(u,Imag(z)); } else { return u + I*v; } }}
/**************************************************************************** double complex clog(double complex z) returns the complex natural logarithm of its argument, using: clog(x + iy) = [ log(x) + 0.5 * log1p(y^2/x^2) ] + I*carg(x + iy) if x > y = [ log(y) + 0.5 * log1p(x^2/y^2) ] + I*carg(x + iy) otherwise the real part is "mathematically" equivalent to log |z|, but the alternative form is used to avoid spurious under/overflow.
Calls: fabs, log1p, log, carg. ****************************************************************************/
double complex clog ( double complex z ) { static const double inf = __builtin_inf(); double large, small, temp; double complex w; long double ratio; Imag(w) = carg(z); // handle x,y = � if ((__builtin_fabs(Real(z)) == inf) || (__builtin_fabs(Imag(z)) == inf)) { Real(w) = inf; return w; } // handle x,y = NaN if (Real(z) != Real(z)) return Real(z) + I*__builtin_copysign(Real(z),Imag(z)); if (Imag(z) != Imag(z)) return Imag(z) + I*Imag(z); large = __builtin_fabs(Real(z)); small = __builtin_fabs(Imag(z)); if (large < small) { temp = large; large = small; small = temp; } Real(w) = log(large); // if small == 0, then Re(clog(z)) = log(large). This sets div-by-zero when appropriate (if large is also 0). if (small == 0.0) return w; // if large == 1 if (large == 1.0) { Real(w) = 0.5*log1p(small*small); // any underflow here is deserved. return w; } // compute small/large in long double to avoid undue underflow. ratio = (long double)small / (long double)large; if (ratio > 0x1.0p-53L) { /* if ratio is smaller than 2^-53, then * 1/2 log1p(ratio^2) ~ 2^-106 < 1/2 an ulp of log(large), so it can't affect the final result. */ Real(w) += 0.5*log1p((double)(ratio*ratio)); } return w;}
float complex clogf ( float complex z ) { static const float inf = __builtin_inff(); float large, small, temp; float complex w; double ratio; Imag(w) = cargf(z); // handle x,y = � if ((__builtin_fabsf(Real(z)) == inf) || (__builtin_fabsf(Imag(z)) == inf)) { Real(w) = inf; return w; } // handle x,y = NaN if (Real(z) != Real(z)) return Real(z) + I*__builtin_copysignf(Real(z),Imag(z)); if (Imag(z) != Imag(z)) return Imag(z) + I*Imag(z); large = __builtin_fabsf(Real(z)); small = __builtin_fabsf(Imag(z)); if (large < small) { temp = large; large = small; small = temp; } Real(w) = logf(large); // if small == 0, then Re(clog(z)) = log(large). This sets div-by-zero when appropriate (if large is also 0). if (small == 0.0f) return w; // if large == 1 if (large == 1.0f) { Real(w) = 0.5f*log1pf(small*small); // underflow here is deserved. return w; } // compute small/large in double to avoid undue underflow. ratio = (double)small / (double)large; if (ratio > 0x1.0p-24) { Real(w) += 0.5f*log1pf((float)(ratio*ratio)); } return w;}
long double complex clogl ( long double complex z ) { static const long double inf = __builtin_infl(); long double x,y; long double complex w; long double ratio; Imag(w) = cargl(z); // handle x,y = � if ((__builtin_fabsl(Real(z)) == inf) || (__builtin_fabsl(Imag(z)) == inf)) { Real(w) = inf; return w; } // handle x,y = NaN if (Real(z) != Real(z)) return Real(z) + I*__builtin_copysignl(Real(z),Imag(z)); if (Imag(z) != Imag(z)) return Imag(z) + I*Imag(z); x = __builtin_fabsl(Real(z)); y = __builtin_fabsl(Imag(z)); ld_parts *large = (ld_parts*) &x; ld_parts *small = (ld_parts*) &y; if (large->ld < small->ld) { ld_parts *p = large; large = small; small = p; } Real(w) = logl(large->ld); // if small == 0, then Re(clog(z)) = log(large). This sets div-by-zero when appropriate (if large is also 0). if (small->ld == 0.0L) return w; // if large == 1 if (large->ld == 1.0L) { Real(w) = 0.5L*log1pl((small->ld)*(small->ld)); // underflow here is deserved. return w; } if (large->sexp - small->sexp < 64) { // if large and small are of roughly comparable magnitude, then the 0.5 * log1p(small^2/large^2) term is // non-negligable. ratio = small->ld / large->ld; Real(w) += 0.5L*log1pl(ratio*ratio); } return w;}
/**************************************************************************** void cosisin(x, complex *z) returns cos(x) + i sin(x) computed using the x87 fsincos instruction. Implemented in s_cosisin.s Called by: cexp, csin, ccos, csinh, and ccosh. ****************************************************************************/void cosisin(double x, double complex *z);void cosisinf(float x, float complex *z);void cosisinl(long double x, long double complex *z);
/**************************************************************************** double complex csin(double complex z) returns the complex sine of its argument. sin(z) = -i sinh(iz)
Calls: csinh****************************************************************************/
double complex csin ( double complex z ) { double complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = csinh(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
float complex csinf ( float complex z ) { float complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = csinhf(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
long double complex csinl ( long double complex z ) { long double complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = csinhl(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
/**************************************************************************** double complex ccos(double complex z) returns the complex cosine of its argument. cos(z) = cosh(iz)
Calls: ccosh****************************************************************************/
double complex ccos ( double complex z ) { double complex iz; Real(iz) = -Imag(z); Imag(iz) = Real(z); return ccosh(iz);}
float complex ccosf ( float complex z ) { float complex iz; Real(iz) = -Imag(z); Imag(iz) = Real(z); return ccoshf(iz);}
long double complex ccosl ( long double complex z ) { long double complex iz; Real(iz) = -Imag(z); Imag(iz) = Real(z); return ccoshl(iz);}
/**************************************************************************** double complex csinh(double complex z) returns the complex hyperbolic sine of its argument. The algorithm is based upon the identity: sinh(x + i*y) = cos(y)*sinh(x) + i*sin(y)*cosh(x). Signaling of spurious overflows, underflows, and invalids is avoided where possible.
Calls: expm1, cosisin****************************************************************************/
double complex csinh ( double complex z ) { static const double INF = __builtin_inf(); double complex w; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysign(Real(z), Imag(z)); return w; } // At this stage, x � NaN. double absx = __builtin_fabs(Real(z)); double reducedx = absx; cosisin(Imag(z), &w); // set w = cos y + i sin y. Real(w) *= __builtin_copysign(1.0, Real(z)); // w = signof(x) cos y + i sin y // Handle x = �� cases. if ((absx == INF) && ((Imag(z) == INF) || (Imag(z) != Imag(z)) || (Imag(z) == 0.0))) { Real(w) = __builtin_copysign(INF, Real(z)); return w; } // Handle x = 0 cases. if (absx == 0.0) { Real(w) = Real(z); // sign of zero needs to be right. return w; } // Argument reduction, if need be. (x is now a finite non-zero number) if ((reducedx < twiceExpOverflowThresh_d) && (reducedx > expOverflowThreshold_d)) { reducedx -= expOverflowThreshold_d; // argument reduction, this is exact. Real(w) *= expOverflowValue_d; // not exact, but good enough. Imag(w) *= expOverflowValue_d; // ditto. } double exm1 = expm1(reducedx); // any overflow here is deserved. if (absx < 0x1p-27) { // |x|^2 is less than an ulp of 1, so only the leading terms of the series for Real(w) *= absx; // cosh = 1 + .... and sinh = x + .... has any effect on the result. } else if (absx > 19.0) { // if |x| > 19, then e^-x is less than an ulp of e^x, so the smaller term in double halfExpX = 0.5 * (exm1 + 1.0); // cosh = (e^x + e^-x) / 2 has no effect, and similarly for Real(w) *= halfExpX; // sinh = (e^x - e^-x) / 2. // only scale the imag part if non-zero (to prevent NaN in overflow*zero) if (Imag(z) != 0.0) Imag(w) *= halfExpX; } else { // the "normal" case, we need to be careful. double twiceExpX = 2.0 * (exm1 + 1.0); /* we use the following to get cosh(x): * * expm1(x)*expm1(x) 2e^x + e^(2x) - 2e^x + 1 e^x + e^-x * 1 + ------------------- = -------------------------- = ------------ = cosh(x) * 2*(1 + expm1(x)) 2e^x 2 */ Imag(w) *= 1.0 + (exm1*exm1)/twiceExpX; /* we use the following to get sinh(x) (up to sign): * * 1 / expm1(x) \ e^(2x) - e^x + e^x - 1 e^x - e^-x * --- | expm1(x) + ------------- | = ------------------------ = ------------ = sinh(x) * 2 \ 1 + expm1(x) / 2e^x 2 */ Real(w) *= 0.5*exm1 + exm1/twiceExpX; } return w;}
float complex csinhf ( float complex z ) { static const float INFf = __builtin_inff(); static const double INF = __builtin_inf(); float complex w; double complex wd; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0f) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysignf(Real(z), Imag(z)); return w; } // At this stage, x � NaN. double absx = (double)__builtin_fabsf(Real(z)); cosisin((double)Imag(z), &wd); // set w = cos y + i sin y. Real(wd) *= __builtin_copysign(1.0, (double)Real(z)); // w = signof(x) cos y + i sin y // Handle x = �� cases. if ((absx == INF) && ((Imag(z) == INFf) || (Imag(z) != Imag(z)) || (Imag(z) == 0.0f))) { Real(w) = __builtin_copysignf(INFf, Real(z)); Imag(w) = (float)Imag(wd); return w; } // Handle x = 0 cases. if (absx == 0.0) { Real(w) = Real(z); // sign of zero needs to be right. Imag(w) = (float)Imag(wd); return w; } double exm1 = expm1(absx); // any overflow here is deserved. if (absx < 0x1p-27) { // |x|^2 is less than an ulp of 1, so only the leading terms of the series for Real(wd) *= absx; // cosh = 1 + .... and sinh = x + .... has any effect on the result. } else if (absx > 19.0) { // if |x| > 19, then e^-x is less than an ulp of e^x, so the smaller term in double halfExpX = 0.5 * (exm1 + 1.0); // cosh = (e^x + e^-x) / 2 has no effect, and similarly for Real(wd) *= halfExpX; // sinh = (e^x - e^-x) / 2. // only scale the imag part if non-zero (to prevent NaN in overflow*zero) if (Imag(z) != 0.0f) Imag(wd) *= halfExpX; } else { // the "normal" case, we need to be careful. double twiceExpX = 2.0 * (exm1 + 1.0); /* we use the following to get cosh(x): * * expm1(x)*expm1(x) 2e^x + e^(2x) - 2e^x + 1 e^x + e^-x * 1 + ------------------- = -------------------------- = ------------ = cosh(x) * 2*(1 + expm1(x)) 2e^x 2 */ Imag(wd) *= 1.0 + (exm1*exm1)/twiceExpX; /* we use the following to get sinh(x) (up to sign): * * 1 / expm1(x) \ e^(2x) - e^x + e^x - 1 e^x - e^-x * --- | expm1(x) + ------------- | = ------------------------ = ------------ = sinh(x) * 2 \ 1 + expm1(x) / 2e^x 2 */ Real(wd) *= 0.5*exm1 + exm1/twiceExpX; } Real(w) = (float)Real(wd); Imag(w) = (float)Imag(wd); return w;}
long double complex csinhl ( long double complex z ) { static const long double INFl = __builtin_infl(); long double complex w; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0L) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysignl(Real(z), Imag(z)); return w; } // At this stage, x � NaN. long double absx = __builtin_fabsl(Real(z)); long double reducedx = absx; cosisinl(Imag(z), &w); // set w = cos y + i sin y. Real(w) *= __builtin_copysignl(1.0L, Real(z)); // w = signof(x) cos y + i sin y // Handle x = �� cases. if ((absx == INFl) && ((Imag(z) == INFl) || (Imag(z) != Imag(z)) || (Imag(z) == 0.0L))) { Real(w) = __builtin_copysignl(INFl, Real(z)); return w; } // Handle x = 0 cases. if (absx == 0.0L) { Real(w) = Real(z); // sign of zero needs to be right. return w; } // Argument reduction, if need be. (x is now a finite non-zero number) if ((reducedx < twiceExpOverflowThresh_ld) && (reducedx > expOverflowThreshold_ld)) { reducedx -= expOverflowThreshold_ld; // argument reduction, this is exact. Real(w) *= expOverflowValue_ld; // not exact, but good enough. Imag(w) *= expOverflowValue_ld; // ditto. } long double exm1 = expm1l(reducedx); // any overflow here is deserved. if (absx < 0x1p-32L) { // |x|^2 is less than an ulp of 1, so only the leading terms of the series for Real(w) *= absx; // cosh = 1 + .... and sinh = x + .... has any effect on the result. } else if (absx > 23L) { // if |x| > 23, then e^-x is less than an ulp of e^x, so the smaller term in long double halfExpX = 0.5L * (exm1 + 1.0L); // cosh = (e^x + e^-x) / 2 has no effect, and similarly for Real(w) *= halfExpX; // sinh = (e^x - e^-x) / 2. // only scale the imag part if non-zero (to prevent NaN in overflow*zero) if (Imag(z) != 0.0L) Imag(w) *= halfExpX; } else { // the "normal" case, we need to be careful. long double twiceExpX = 2.0L * (exm1 + 1.0L); /* we use the following to get cosh(x): * * expm1(x)*expm1(x) 2e^x + e^(2x) - 2e^x + 1 e^x + e^-x * 1 + ------------------- = -------------------------- = ------------ = cosh(x) * 2*(1 + expm1(x)) 2e^x 2 */ Imag(w) *= 1.0L + (exm1*exm1)/twiceExpX; /* we use the following to get sinh(x) (up to sign): * * 1 / expm1(x) \ e^(2x) - e^x + e^x - 1 e^x - e^-x * --- | expm1(x) + ------------- | = ------------------------ = ------------ = sinh(x) * 2 \ 1 + expm1(x) / 2e^x 2 */ Real(w) *= 0.5L*exm1 + exm1/twiceExpX; } return w;}
/**************************************************************************** double complex ccosh(double complex z) returns the complex hyperbolic cosine of its argument. The algorithm is based upon the identity: cosh(x + i*y) = cos(y)*cosh(x) + i*sin(y)*sinh(x). Signaling of spurious overflows, underflows, and invalids is avoided where possible.
Calls: expm1, cosisin****************************************************************************/
double complex ccosh ( double complex z ) { static const double INF = __builtin_inf(); double complex w; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysign(Real(z), Imag(z)); return w; } // At this stage, x � NaN. double absx = __builtin_fabs(Real(z)); double reducedx = absx; cosisin(Imag(z), &w); // set w = cos y + i sin y. Imag(w) *= __builtin_copysign(1.0, Real(z)); // w = cos y + i sin y * signof(x) // Handle x = �� cases. if ((absx == INF) && ((Imag(z) == INF) || (Imag(z) != Imag(z)) || (Imag(z) == 0.0))) { Real(w) = INF; return w; } // Handle x = 0 cases. if (absx == 0.0) { Imag(w) = Real(z) * __builtin_copysign(1.0, Imag(z)); // finesse the sign of zero. return w; } // Argument reduction, if need be. (x is now a finite non-zero number) if ((reducedx < twiceExpOverflowThresh_d) && (reducedx > expOverflowThreshold_d)) { reducedx -= expOverflowThreshold_d; // argument reduction, this is exact. Real(w) *= expOverflowValue_d; // not exact, but good enough. Imag(w) *= expOverflowValue_d; // ditto. } double exm1 = expm1(reducedx); // any overflow here is deserved. if (absx < 0x1p-27) { // |x|^2 is less than an ulp of 1, so only the leading terms of the series for Imag(w) *= absx; // cosh = 1 + .... and sinh = x + .... has any effect on the result. } else if (absx > 19.0) { // if |x| > 19, then e^-x is less than an ulp of e^x, so the smaller term in double halfExpX = 0.5 * (exm1 + 1.0); // cosh = (e^x + e^-x) / 2 has no effect, and similarly for Real(w) *= halfExpX; // sinh = (e^x - e^-x) / 2. // only scale the imag part if non-zero (to prevent NaN in overflow*zero) if (Imag(z) != 0.0) Imag(w) *= halfExpX; } else { // the "normal" case, we need to be careful. double twiceExpX = 2.0 * (exm1 + 1.0); /* we use the following to get cosh(x): * * expm1(x)*expm1(x) 2e^x + e^(2x) - 2e^x + 1 e^x + e^-x * 1 + ------------------- = -------------------------- = ------------ = cosh(x) * 2*(1 + expm1(x)) 2e^x 2 */ Real(w) *= 1.0 + (exm1*exm1)/twiceExpX; /* we use the following to get sinh(x) (up to sign): * * 1 / expm1(x) \ e^(2x) - e^x + e^x - 1 e^x - e^-x * --- | expm1(x) + ------------- | = ------------------------ = ------------ = sinh(x) * 2 \ 1 + expm1(x) / 2e^x 2 */ Imag(w) *= 0.5*exm1 + exm1/twiceExpX; } return w;}
float complex ccoshf ( float complex z ) { static const float INFf = __builtin_inff(); static const double INF = __builtin_inf(); double complex wd; float complex w; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0f) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysignf(Real(z), Imag(z)); return w; } // At this stage, x � NaN. double absx = (double)__builtin_fabsf(Real(z)); cosisin((double)Imag(z), &wd); // set w = cos y + i sin y. Imag(wd) *= __builtin_copysign(1.0, (double)Real(z)); // w = cos y + i sin y * signof(x) // Handle x = �� cases. if ((absx == INF) && ((Imag(z) == INFf) || (Imag(z) != Imag(z)) || (Imag(z) == 0.0f))) { Real(w) = INFf; Imag(w) = (float)Imag(wd); return w; } // Handle x = 0 cases. if (absx == 0.0) { Imag(w) = Real(z) * __builtin_copysignf(1.0f, Imag(z)); // finesse the sign of zero. Real(w) = (float)Real(wd); return w; } double exm1 = expm1(absx); // any overflow here is deserved. if (absx < 0x1p-27) { // |x|^2 is less than an ulp of 1, so only the leading terms of the series for Imag(wd) *= absx; // cosh = 1 + .... and sinh = x + .... has any effect on the result. } else if (absx > 19.0) { // if |x| > 19, then e^-x is less than an ulp of e^x, so the smaller term in double halfExpX = 0.5 * (exm1 + 1.0); // cosh = (e^x + e^-x) / 2 has no effect, and similarly for Real(wd) *= halfExpX; // sinh = (e^x - e^-x) / 2. // only scale the imag part if non-zero (to prevent NaN in overflow*zero) if (Imag(z) != 0.0) Imag(wd) *= halfExpX; } else { // the "normal" case, we need to be careful. double twiceExpX = 2.0 * (exm1 + 1.0); /* we use the following to get cosh(x): * * expm1(x)*expm1(x) 2e^x + e^(2x) - 2e^x + 1 e^x + e^-x * 1 + ------------------- = -------------------------- = ------------ = cosh(x) * 2*(1 + expm1(x)) 2e^x 2 */ Real(wd) *= 1.0 + (exm1*exm1)/twiceExpX; /* we use the following to get sinh(x) (up to sign): * * 1 / expm1(x) \ e^(2x) - e^x + e^x - 1 e^x - e^-x * --- | expm1(x) + ------------- | = ------------------------ = ------------ = sinh(x) * 2 \ 1 + expm1(x) / 2e^x 2 */ Imag(wd) *= 0.5*exm1 + exm1/twiceExpX; } Real(w) = (float)Real(wd); Imag(w) = (float)Imag(wd); return w;}
long double complex ccoshl ( long double complex z ) { static const long double INFl = __builtin_infl(); long double complex w; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0L) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysignl(Real(z), Imag(z)); return w; } // At this stage, x � NaN. long double absx = __builtin_fabsl(Real(z)); long double reducedx = absx; cosisinl(Imag(z), &w); // set w = cos y + i sin y. Imag(w) *= __builtin_copysignl(1.0, Real(z)); // w = cos y + i sin y * signof(x) // Handle x = �� cases. if ((absx == INFl) && ((Imag(z) == INFl) || (Imag(z) != Imag(z)) || (Imag(z) == 0.0L))) { Real(w) = INFl; return w; } // Handle x = 0 cases. if (absx == 0.0L) { Imag(w) = Real(z) * __builtin_copysignl(1.0, Imag(z)); // finesse the sign of zero. return w; } // Argument reduction, if need be. (x is now a finite non-zero number) if ((reducedx < twiceExpOverflowThresh_ld) && (reducedx > expOverflowThreshold_ld)) { reducedx -= expOverflowThreshold_ld; // argument reduction, this is exact. Real(w) *= expOverflowValue_ld; // not exact, but good enough. Imag(w) *= expOverflowValue_ld; // ditto. } long double exm1 = expm1l(reducedx); // any overflow here is deserved. if (absx < 0x1p-32L) { // |x|^2 is less than an ulp of 1, so only the leading terms of the series for Imag(w) *= absx; // cosh = 1 + .... and sinh = x + .... has any effect on the result. } else if (absx > 23L) { // if |x| > 23, then e^-x is less than an ulp of e^x, so the smaller term in long double halfExpX = 0.5L * (exm1 + 1.0L); // cosh = (e^x + e^-x) / 2 has no effect, and similarly for Real(w) *= halfExpX; // sinh = (e^x - e^-x) / 2. // only scale the imag part if non-zero (to prevent NaN in overflow*zero) if (Imag(z) != 0.0L) Imag(w) *= halfExpX; } else { // the "normal" case, we need to be careful. long double twiceExpX = 2.0L * (exm1 + 1.0L); /* we use the following to get cosh(x): * * expm1(x)*expm1(x) 2e^x + e^(2x) - 2e^x + 1 e^x + e^-x * 1 + ------------------- = -------------------------- = ------------ = cosh(x) * 2*(1 + expm1(x)) 2e^x 2 */ Real(w) *= 1.0L + (exm1*exm1)/twiceExpX; /* we use the following to get sinh(x) (up to sign): * * 1 / expm1(x) \ e^(2x) - e^x + e^x - 1 e^x - e^-x * --- | expm1(x) + ------------- | = ------------------------ = ------------ = sinh(x) * 2 \ 1 + expm1(x) / 2e^x 2 */ Imag(w) *= 0.5L*exm1 + exm1/twiceExpX; } return w;}
/**************************************************************************** double complex cexp(double complex z) returns the complex exponential of its argument. The algorithm is based upon the identity: exp(x + i*y) = cos(y)*exp(x) + i*sin(y)*exp(x). Signaling of spurious overflows and invalids is avoided where possible. CONSTANTS: expOverflowValue_d = exp(709.0) to double precision
Calls: cosisin and exp.****************************************************************************/
double complex cexp ( double complex z ) { static const double INF = __builtin_inf(); double complex w; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysign(Real(z), Imag(z)); return w; } // Handle x = -�, y = � or NaN: if ((Real(z) == -INF) && ((__builtin_fabs(Imag(z)) == INF) || (Imag(z) != Imag(z)))) { Real(w) = 0.0; Imag(w) = __builtin_copysign(0.0, Imag(z)); return w; } if (Imag(z) == 0.0) { // exact exp(x + 0i) case. Real(w) = exp(Real(z)); Imag(w) = __builtin_copysign(0.0, Imag(z)); return w; } // At this stage, x � NaN, and extraordinary x = -� cases are sorted. y � 0. cosisin(Imag(z), &w); // set w = cos(y) + i sin(y) // Handle x = +� cases. if ((Real(z) == INF) && ((Imag(z) == INF) || (Imag(z) != Imag(z)))) { Real(w) = INF; // cexp(� + yi) = � + NaNi for y = NaN or �. return w; // cexp(� + yi) for finite y falls through. } // At this point, x � NaN, +inf, y � 0, and all remaining cases fall through double x = Real(z); if ((x < twiceExpOverflowThresh_d) && (x > expOverflowThreshold_d)) { x -= expOverflowThreshold_d; // argument reduction, this is exact. Real(w) *= expOverflowValue_d; // not exact, but good enough. Imag(w) *= expOverflowValue_d; // ditto. } double scale = exp(x); Real(w) *= scale; Imag(w) *= scale; return w;}
float complex cexpf ( float complex z ) { static const float INFf = __builtin_inff(); float complex w; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0f) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysignf(Real(z), Imag(z)); return w; } // Handle x = -�, y = � or NaN: if ((Real(z) == -INFf) && ((__builtin_fabsf(Imag(z)) == INFf) || (Imag(z) != Imag(z)))) { Real(w) = 0.0f; Imag(w) = __builtin_copysignf(0.0f, Imag(z)); return w; } if (Imag(z) == 0.0f) { // exact exp(x + 0i) case. Real(w) = expf(Real(z)); Imag(w) = __builtin_copysignf(0.0f, Imag(z)); return w; } double complex wd; // At this stage, x � NaN, and extraordinary x = -� cases are sorted. y � 0. cosisin((double)Imag(z), &wd); // set w = cos(y) + i sin(y) // Handle x = +� cases. if ((Real(z) == INFf) && ((Imag(z) == INFf) || (Imag(z) != Imag(z)))) { Real(w) = INFf; // cexp(� + yi) = � + NaNi for y = NaN or �. Imag(w) = (float)Imag(wd); return w; // cexp(� + yi) for finite y falls through. } // At this point, x � NaN, +inf, y � 0, and all remaining cases fall through
double scale = exp((double)Real(z)); Real(w) = (float)(scale*Real(wd)); Imag(w) = (float)(scale*Imag(wd)); return w;}
long double complex cexpl ( long double complex z ) { static const long double INFl = __builtin_infl(); long double complex w; // Handle x = NaN first if (Real(z) != Real(z)) { Real(w) = Real(z); if (Imag(z) == 0.0L) // cexp(NaN + 0i) = NaN + 0i Imag(w) = Imag(z); else // cexp(NaN + yi) = NaN + NaNi, for y � 0 Imag(w) = __builtin_copysignl(Real(z), Imag(z)); return w; } // Handle x = -�, y = � or NaN: if ((Real(z) == -INFl) && ((__builtin_fabsl(Imag(z)) == INFl) || (Imag(z) != Imag(z)))) { Real(w) = 0.0L; Imag(w) = __builtin_copysignl(0.0L, Imag(z)); return w; } if (Imag(z) == 0.0L) { // exact exp(x + 0i) case. Real(w) = expl(Real(z)); Imag(w) = __builtin_copysignl(0.0L, Imag(z)); return w; } // At this stage, x � NaN, and extraordinary x = -� cases are sorted. y � 0. cosisinl(Imag(z), &w); // set w = cos(y) + i sin(y) // Handle x = +� cases. if ((Real(z) == INFl) && ((Imag(z) == INFl) || (Imag(z) != Imag(z)))) { Real(w) = INFl; // cexp(� + yi) = � + NaNi for y = NaN or �, � + 0i for y = 0. return w; // cexp(� + yi) for finite y falls through. } // At this point, x � NaN, +inf, y � 0, and all remaining cases fall through long double x = Real(z); if ((x < twiceExpOverflowThresh_ld) && (x > expOverflowThreshold_ld)) { x -= expOverflowThreshold_ld; // argument reduction, this is exact. Real(w) *= expOverflowValue_ld; // not exact, but good enough. Imag(w) *= expOverflowValue_ld; // ditto. } long double scale = expl(x); Real(w) *= scale; Imag(w) *= scale; return w;} /**************************************************************************** double complex cpow(double complex x, double complex y) returns the complex result of x^y. The algorithm is based upon the identity: x^y = cexp(y*clog(x)). Calls: clog, cexp.****************************************************************************/
double complex cpow ( double complex x, double complex y ) /* (complex x)^(complex y) */{ double complex logval,z; logval = clog(x); /* complex logarithm of x */ Real(z) = Real(y)*Real(logval) - Imag(y)*Imag(logval); /* multiply by y */ Imag(z) = Real(y)*Imag(logval) + Imag(y)*Real(logval); return (cexp(z)); /* return complex exponential */}
float complex cpowf ( float complex x, float complex y ) /* (complex x)^(complex y) */{ float complex logval,z; logval = clogf(x); /* complex logarithm of x */ Real(z) = Real(y)*Real(logval) - Imag(y)*Imag(logval); /* multiply by y */ Imag(z) = Real(y)*Imag(logval) + Imag(y)*Real(logval); return (cexpf(z)); /* return complex exponential */}
long double complex cpowl ( long double complex x, long double complex y ) /* (complex x)^(complex y) */{ long double complex logval,z; logval = clogl(x); /* complex logarithm of x */ Real(z) = Real(y)*Real(logval) - Imag(y)*Imag(logval); /* multiply by y */ Imag(z) = Real(y)*Imag(logval) + Imag(y)*Real(logval); return (cexpl(z)); /* return complex exponential */}
/**************************************************************************** double complex ctanh(double complex x) returns the complex hyperbolic tangent of its argument. The algorithm is from Kahan's paper and is based on the identity: tanh(x+i*y) = (sinh(2*x) + i*sin(2*y))/(cosh(2*x) + cos(2*y)) = (cosh(x)*sinh(x)*cscsq + i*tan(y))/(1+cscsq*sinh(x)*sinh(x)), where cscsq = 1/(cos(y)*cos(y). For large values of ze.re, spurious overflow and invalid signaling is avoided. CONSTANTS: FPKASINHOM4 = asinh(nextafterd(+infinity,0.0))/4.0 to double precision FPKINF = +infinity Calls: tan, sinh, sqrt, fabs, feholdexcept, feraiseexcept, feclearexcept, and feupdateenv.****************************************************************************/
double complex ctanh( double complex z ){ static const double INF = __builtin_inf(); double x = __builtin_fabs(Real(z)); double y = __builtin_fabs(Imag(z)); double sinhval, coshval, tanval, exm1, cscsq; double complex w; if (x == INF) { w = 1.0 + I*__builtin_copysign(0.0, sin(2.0*y)); // ctanh(� + iy) = 1.0 � i0 } else if (Imag(z) != Imag(z) || Real(z) != Real(z)) { if (Imag(z) == 0.0) { w = Real(z) + I*0.0; // ctanh(NaN + i0) = NaN + i0 } else { Real(w) = Real(z) + Imag(z); // ctanh(NaN) = NaN + iNaN Imag(w) = Real(w); } } else if (y == INF) { Real(w) = y - y; // ctanh(x + i�) = NaN + iNaN (invalid) Imag(w) = Real(w); } else if (x > 19.0) { w = 1.0 + I*__builtin_copysign(0.0, sin(2.0*y)); // if x is big, tanh(z) = 1 � i0 } else { // edge case free! tanval = tan(y); cscsq = 1.0 + tanval*tanval; // cscsq = 1/cos^2(y) if (x < 0x1p-27) { coshval = 1.0; sinhval = x; } else { exm1 = expm1(x); coshval = 1.0 + 0.5*(exm1*exm1)/(exm1 + 1.0); sinhval = 0.5*(exm1 + exm1/(exm1 + 1.0)); } Real(w) = cscsq * coshval * sinhval / (1.0 + cscsq * sinhval * sinhval); Imag(w) = tanval / (1.0 + cscsq * sinhval * sinhval); } // Patch up signs of return value Real(w) = __builtin_copysign(Real(w),Real(z)); Imag(w) *= __builtin_copysign(1.0,Imag(z)); return w;}
float complex ctanhf( float complex z ){ static const float INFf = __builtin_inff(); float x = __builtin_fabsf(Real(z)); float y = __builtin_fabsf(Imag(z)); double sinhval, coshval, tanval, exm1, cscsq; float complex w; if (x == INFf) { w = 1.0f + I*__builtin_copysignf(0.0f, sinf(2.0f*y)); // ctanh(� + iy) = 1.0 � i0 } else if (Imag(z) != Imag(z) || Real(z) != Real(z)) { if (Imag(z) == 0.0f) { w = Real(z) + I*0.0f; // ctanh(NaN + i0) = NaN + i0 } else { Real(w) = Real(z) + Imag(z); // ctanh(NaN) = NaN + iNaN Imag(w) = Real(w); } } else if (y == INFf) { Real(w) = y - y; // ctanh(x + i�) = NaN + iNaN (invalid) Imag(w) = Real(w); } else if (x > 19.0f) { w = 1.0f + I*__builtin_copysignf(0.0f, sinf(2.0f*y)); // if x is big, tanh(z) = 1 � i0 } else { // edge case free! tanval = (double)tanf(y); cscsq = 1.0 + tanval*tanval; // cscsq = 1/cos^2(y) if (x < 0x1p-13f) { coshval = 1.0; sinhval = x; } else { exm1 = (double)expm1f(x); coshval = 1.0 + 0.5*(exm1*exm1)/(exm1 + 1.0); sinhval = 0.5*(exm1 + exm1/(exm1 + 1.0)); } Real(w) = (float)(cscsq * coshval * sinhval / (1.0 + cscsq * sinhval * sinhval)); Imag(w) = (float)(tanval / (1.0 + cscsq * sinhval * sinhval)); } // Patch up signs of return value Real(w) = __builtin_copysignf(Real(w),Real(z)); Imag(w) *= __builtin_copysignf(1.0f,Imag(z)); return w;}
long double complex ctanhl( long double complex z ){ static const long double INFl = __builtin_infl(); long double x = __builtin_fabsl(Real(z)); long double y = __builtin_fabsl(Imag(z)); long double sinhval, coshval, tanval, exm1, cscsq; long double complex w; if (x == INFl) { w = 1.0l + I*__builtin_copysignl(0.0l, sinl(2.0l*y)); // ctanh(� + iy) = 1.0 � i0 } else if (Imag(z) != Imag(z) || Real(z) != Real(z)) { if (Imag(z) == 0.0l) { w = Real(z) + I*0.0l; // ctanh(NaN + i0) = NaN + i0 } else { Real(w) = Real(z) + Imag(z); // ctanh(NaN) = NaN + iNaN Imag(w) = Real(w); } } else if (y == INFl) { Real(w) = y - y; // ctanh(x + i�) = NaN + iNaN (invalid) Imag(w) = Real(w); } else if (x > 22.0l) { w = 1.0l + I*__builtin_copysignl(0.0l, sinl(2.0l*y)); // if x is big, tanh(z) = 1 � i0 } else { // edge case free! tanval = tanl(y); cscsq = 1.0l + tanval*tanval; // cscsq = 1/cos^2(y) if (x < 0x1p-32l) { coshval = 1.0l; sinhval = x; } else { exm1 = expm1l(x); coshval = 1.0l + 0.5l*(exm1*exm1)/(exm1 + 1.0l); sinhval = 0.5l*(exm1 + exm1/(exm1 + 1.0l)); } Real(w) = cscsq * coshval * sinhval / (1.0l + cscsq * sinhval * sinhval); Imag(w) = tanval / (1.0l + cscsq * sinhval * sinhval); } // Patch up signs of return value Real(w) = __builtin_copysignl(Real(w),Real(z)); Imag(w) *= __builtin_copysignl(1.0l,Imag(z)); return w;}
/**************************************************************************** double complex ctan(double complex x) returns the complex hyperbolic tangent of its argument. Per C99, i ctan(z) = ctanh(iz) Calls: ctanh****************************************************************************/
double complex ctan( double complex z ){ double complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = ctanh(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
float complex ctanf( float complex z ){ float complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = ctanhf(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
long double complex ctanl( long double complex z ){ long double complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = ctanhl(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;} /**************************************************************************** double complex casin(double complex z) returns the complex inverse sine of its argument. The algorithm is from Kahan's paper and is based on the formulae: real(casin(z)) = atan (real(z)/real(csqrt(1.0-z)*csqrt(1.0+z))) imag(casin(z)) = arcsinh(imag(csqrt(1.0-cconj(z))*csqrt(1.0+z))) Calls: arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv.****************************************************************************/
double complex casin ( double complex z ){ double complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = casinh(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
float complex casinf ( float complex z ){ float complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = casinhf(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
long double complex casinl ( long double complex z ){ long double complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = casinhl(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
/**************************************************************************** double complex casinh(double complex z) returns the complex inverse hyperbolic sine of its argument. We compute the function only in the upper-right quadrant of the complex plane, and use the facts that casinh(conj(z)) = conj(casinh(z)) and casinh is odd to get the values on the rest of the plane. within the upper-right quadrant, we use:
casinh(z) = z if |z| is small, ln(2z) if |z| is big, and a rather complicated expression for other values of z. Calls: asinh, csqrt, atan2.****************************************************************************/
double complex casinh ( double complex z ) { static const double INF = __builtin_inf(); static const double ln2 = 0x1.62e42fefa39ef358p-1; static const double sqrt1_2 = 0x1.6a09e667f3bcc908p-1; double complex w; double x = __builtin_fabs(Real(z)); double y = __builtin_fabs(Imag(z)); double u, xSquared, tmp; double complex sqrt1Plusiz, sqrt1PlusizBar; // If |z| == inf, then casinh(z) = inf + carg(z) if ((x == INF) || (y == INF)) { Real(w) = INF; Imag(w) = atan2(y,x); } // If z = NaN, casinh(z) = NaN, with the special case that casinh(NaN + i0) = NaN + i0. else if ((x != x) || (y != y)) { if (y == 0.0) w = z; else { Real(w) = x + y; Imag(w) = x + y; } }
// at this point x,y are finite, non-nan. else { // If z is small, then casinh(z) = z - z^3/6 + ... = z within half an ulp if ((x < 0x1p-27) && (y < 0x1p-27)) { Real(w) = x; Imag(w) = y; } // If z is big, then casinh(z) = log2 + log(z) + terms smaller than half an ulp else if ((x > 0x1p27) || (y > 0x1p27)) { w = clog(x + I*y); Real(w) += ln2; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. */ else { // Compute sqrt1Plusiz = sqrt(1-y + ix) u = 1.0 - y; xSquared = (x < 0x1p-106 ? 0.0 : x*x); // Avoid underflows. Faster via mask? if (u == 0.0) { Real(sqrt1Plusiz) = sqrt1_2 * __builtin_sqrt(x); // Avoid spurious underflows in this case Imag(sqrt1Plusiz) = Real(sqrt1Plusiz); // by using the simpler formula. } else { // No underflow or overflow is possible. Real(sqrt1Plusiz) = __builtin_sqrt(0.5*(__builtin_sqrt(u*u + xSquared) + __builtin_fabs(u))); tmp = 0.5 * (x / Real(sqrt1Plusiz)); if (u < 0.0) { Imag(sqrt1Plusiz) = Real(sqrt1Plusiz); Real(sqrt1Plusiz) = tmp; } else { Imag(sqrt1Plusiz) = tmp; } } // Compute sqrt1PlusizBar = sqrt(1+y + ix). No underflow or overflow is possible. u = 1.0 + y; Real(sqrt1PlusizBar) = __builtin_sqrt(0.5*(__builtin_sqrt(u*u + xSquared) + u)); Imag(sqrt1PlusizBar) = x / (2.0*Real(sqrt1PlusizBar)); // Magic formulas from Kahan. Real(w) = asinh(Real(sqrt1Plusiz)*Imag(sqrt1PlusizBar) + Imag(sqrt1Plusiz)*Real(sqrt1PlusizBar)); Imag(w) = atan2(y, Real(sqrt1Plusiz)*Real(sqrt1PlusizBar) + Imag(sqrt1Plusiz)*Imag(sqrt1PlusizBar)); } } // Patch up signs to handle z in quadrants II - IV, using symmetry. Real(w) = __builtin_copysign(Real(w), Real(z)); Imag(w) = __builtin_copysign(Imag(w), Imag(z)); return w;}
float complex casinhf ( float complex z ) { static const float INFf = __builtin_inff(); static const float ln2f = 0x1.62e42fefa39ef358p-1f; static const float sqrt1_2f = 0x1.6a09e667f3bcc908p-1f; float complex w; float x = __builtin_fabsf(Real(z)); float y = __builtin_fabsf(Imag(z)); float u, xSquared, tmp; float complex sqrt1Plusiz, sqrt1PlusizBar; // If |z| == inf, then casinh(z) = inf + carg(z) if ((x == INFf) || (y == INFf)) { Real(w) = INFf; Imag(w) = atan2f(y,x); } // If z = NaN, casinh(z) = NaN, with the special case that casinh(NaN + i0) = NaN + i0. else if ((x != x) || (y != y)) { if (y == 0.0f) w = z; else { Real(w) = x + y; Imag(w) = x + y; } } // at this point x,y are finite, non-nan. else { // If z is small, then casinhf(z) = z - z^3/6 + ... = z within half an ulp if ((x < 0x1p-13f) && (y < 0x1p-13f)) { Real(w) = x; Imag(w) = y; } // If z is big, then casinh(z) = log2 + log(z) + terms smaller than half an ulp else if ((x > 0x1p13f) || (y > 0x1p13f)) { w = clogf(x + I*y); Real(w) += ln2f; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. */ else { // Compute sqrt1Plusiz = sqrt(1-y + ix) u = 1.0f - y; xSquared = (x < 0x1p-52f ? 0.0f : x*x); // Avoid underflows. Faster via mask? if (u == 0.0f) { Real(sqrt1Plusiz) = sqrt1_2f * __builtin_sqrtf(x); // Avoid spurious underflows in this case Imag(sqrt1Plusiz) = Real(sqrt1Plusiz); // by using the simpler formula. } else { // No underflow or overflow is possible. Real(sqrt1Plusiz) = __builtin_sqrtf(0.5f*(__builtin_sqrtf(u*u + xSquared) + __builtin_fabsf(u))); tmp = 0.5f * (x / Real(sqrt1Plusiz)); if (u < 0.0f) { Imag(sqrt1Plusiz) = Real(sqrt1Plusiz); Real(sqrt1Plusiz) = tmp; } else { Imag(sqrt1Plusiz) = tmp; } } // Compute sqrt1PlusizBar = sqrt(1+y + ix). No underflow or overflow is possible. u = 1.0f + y; Real(sqrt1PlusizBar) = __builtin_sqrtf(0.5f*(__builtin_sqrtf(u*u + xSquared) + u)); Imag(sqrt1PlusizBar) = x / (2.0f*Real(sqrt1PlusizBar)); // Magic formulas from Kahan. Real(w) = asinhf(Real(sqrt1Plusiz)*Imag(sqrt1PlusizBar) + Imag(sqrt1Plusiz)*Real(sqrt1PlusizBar)); Imag(w) = atan2f(y, Real(sqrt1Plusiz)*Real(sqrt1PlusizBar) + Imag(sqrt1Plusiz)*Imag(sqrt1PlusizBar)); } } // Patch up signs to handle z in quadrants II - IV, using symmetry. Real(w) = __builtin_copysignf(Real(w), Real(z)); Imag(w) = __builtin_copysignf(Imag(w), Imag(z)); return w;}
long double complex casinhl ( long double complex z ) { static const long double INFl = __builtin_infl(); static const long double ln2l = 0x1.62e42fefa39ef358p-1L; static const long double sqrt1_2l = 0x1.6a09e667f3bcc908p-1L; long double complex w; long double x = __builtin_fabsl(Real(z)); long double y = __builtin_fabsl(Imag(z)); long double u, xSquared, tmp; long double complex sqrt1Plusiz, sqrt1PlusizBar; // If |z| == inf, then casinh(z) = inf + carg(z) if ((x == INFl) || (y == INFl)) { Real(w) = INFl; Imag(w) = atan2l(y,x); } // If z = NaN, casinh(z) = NaN, with the special case that casinh(NaN + i0) = NaN + i0. else if ((x != x) || (y != y)) { if (y == 0.0l) w = z; else { Real(w) = x + y; Imag(w) = x + y; } } // at this point x,y are finite, non-nan. else { // If z is small, then casinhl(z) = z - z^3/6 + ... = z within half an ulp if ((x < 0x1p-32l) && (y < 0x1p-32l)) { Real(w) = x; Imag(w) = y; } // If z is big, then casinhl(z) = log2 + log(z) + terms smaller than half an ulp else if ((x > 0x1p32l) || (y > 0x1p32l)) { w = clogl(x + I*y); Real(w) += ln2l; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. */ else { u = 1.0l - y; xSquared = (x < 0x1p-128l ? 0.0l : x*x); // Avoid underflows. Faster via mask? if (u == 0.0l) { Real(sqrt1Plusiz) = sqrt1_2l * __builtin_sqrtl(x); // Avoid spurious underflows in this case Imag(sqrt1Plusiz) = Real(sqrt1Plusiz); // by using the simpler formula. } else { // No underflow or overflow is possible. Real(sqrt1Plusiz) = __builtin_sqrtl(0.5l*(__builtin_sqrtl(u*u + xSquared) + __builtin_fabsl(u))); tmp = 0.5 * (x / Real(sqrt1Plusiz)); if (u < 0.0l) { Imag(sqrt1Plusiz) = Real(sqrt1Plusiz); Real(sqrt1Plusiz) = tmp; } else { Imag(sqrt1Plusiz) = tmp; } } // Compute sqrt1PlusizBar = sqrt(1+y + ix). No underflow or overflow is possible. u = 1.0l + y; Real(sqrt1PlusizBar) = __builtin_sqrtl(0.5l*(__builtin_sqrtl(u*u + xSquared) + u)); Imag(sqrt1PlusizBar) = x / (2.0l*Real(sqrt1PlusizBar)); // Magic formulas from Kahan. Real(w) = asinhl(Real(sqrt1Plusiz)*Imag(sqrt1PlusizBar) + Imag(sqrt1Plusiz)*Real(sqrt1PlusizBar)); Imag(w) = atan2l(y, Real(sqrt1Plusiz)*Real(sqrt1PlusizBar) + Imag(sqrt1Plusiz)*Imag(sqrt1PlusizBar)); } } // Patch up signs to handle z in quadrants II - IV, using symmetry. Real(w) = __builtin_copysignl(Real(w), Real(z)); Imag(w) = __builtin_copysignl(Imag(w), Imag(z)); return w;}
/**************************************************************************** double complex cacos(double complex z) returns the complex inverse cosine of its argument. The algorithm is from Kahan's paper and is based on the formulae: real(cacos(z)) = 2.0*atan(real(csqrt(1.0-z)/real(csqrt(1.0+z)))) imag(cacos(z)) = arcsinh(imag(csqrt(1.0-z)*csqrt(cconj(1.0+z)))) Calls: arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv.****************************************************************************/
double complex cacos ( double complex z ){ static const double INF = __builtin_inf(); static const double ln2 = 0x1.62e42fefa39ef358p-1; static const double sqrt1_2 = 0x1.6a09e667f3bcc908p-1; static const double pi2 = 0x1.921fb54442d1846ap0;
double complex w; double x = __builtin_fabs(Real(z)); double y = __builtin_fabs(Imag(z)); double u, ySquared, tmp; double complex sqrt1Plusz, sqrt1Minusz; // If |z| == inf, then cacos(z) = carg(z) - inf * I if ((x == INF) || (y == INF)) { Imag(w) = -INF; Real(w) = atan2(y,x); } // If z = NaN, cacos(z) = NaN, with the special case that cacos(0 + iNaN) = �/2 + iNaN. else if ((x != x) || (y != y)) { if (x == 0.0) Real(w) = pi2; else Real(w) = x + y; Imag(w) = x + y; } // at this point x,y are finite, non-nan. else { // If z is small, then cacos(z) = �/2 - z + z^3/6 + ... = �/2 - z within half an ulp if ((x < 0x1p-27) && (y < 0x1p-27)) { Real(w) = pi2 - x; Imag(w) = -y; } // If z is big, then cacos(z) = -i * (log2 + log(z)) + terms smaller than half an ulp else if ((x > 0x1p27) || (y > 0x1p27)) { w = clog(x + I*y) + ln2; const double tmp = __real__ w; __real__ w = __imag__ w; __imag__ w = -tmp; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Real(w) = 2*atan2( Re(sqrt(1-z)), Re(sqrt(1+z)) ) * Imag(w) = asinh( Im( sqrt(1-z)*sqrt(1+conj(z)) ) ) * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. Analysis of roundoff is a bit harder, but goes though just fine. */ else { ySquared = (y < 0x1p-106 ? 0.0 : y*y); // Avoid underflows. Faster via mask? // Compute sqrt1Plusz = sqrt(1+x + iy) u = 1.0 + x; Real(sqrt1Plusz) = __builtin_sqrt(0.5*(__builtin_sqrt(u*u + ySquared) + u)); Imag(sqrt1Plusz) = 0.5 * (y / Real(sqrt1Plusz)); // Compute sqrt1Minusz = sqrt(1-x - iy) u = 1.0 - x; if (u == 0.0) { Real(sqrt1Minusz) = sqrt1_2 * __builtin_sqrt(y); // Avoid spurious underflows in this case Imag(sqrt1Minusz) = -Real(sqrt1Minusz); // by using the simpler formula. } else { // No underflow or overflow is possible. Real(sqrt1Minusz) = __builtin_sqrt(0.5*(__builtin_sqrt(u*u + ySquared) + __builtin_fabs(u))); tmp = 0.5 * (y / Real(sqrt1Minusz)); if (u < 0.0) { Imag(sqrt1Minusz) = -Real(sqrt1Minusz); Real(sqrt1Minusz) = tmp; } else { Imag(sqrt1Minusz) = -tmp; } } // Magic formulas from Kahan. Real(w) = 2.0 * atan2(Real(sqrt1Minusz), Real(sqrt1Plusz)); Imag(w) = asinh( Real(sqrt1Plusz)*Imag(sqrt1Minusz) - Imag(sqrt1Plusz)*Real(sqrt1Minusz) ); } } // Patch up signs to handle z in quadrants II, III & IV, using symmetry. Imag(w) = __builtin_copysign(Imag(w), -Imag(z)); if (Real(z) < 0.0) Real(w) = 2.0 * pi2 - Real(w); // No undue cancellation is possible here - Real(w) < �/2. return w;} float complex cacosf ( float complex z ){ static const float INFf = __builtin_inff(); static const float ln2f = 0x1.62e42fefa39ef358p-1f; static const float pi2f = 0x1.921fb54442d1846ap0f; static const float sqrt1_2f = 0x1.6a09e667f3bcc908p-1f; float complex w; float x = __builtin_fabsf(Real(z)); float y = __builtin_fabsf(Imag(z)); float u, ySquared, tmp; float complex sqrt1Plusz, sqrt1Minusz; // If |z| == inf, then cacos(z) = carg(z) - inf i if ((x == INFf) || (y == INFf)) { Imag(w) = -INFf; Real(w) = atan2f(y,x); } // If z = NaN, cacos(z) = NaN, with the special case that cacos(0 + iNaN) = �/2 + iNaN. else if ((x != x) || (y != y)) { if (x == 0.0f) Real(w) = pi2f; else Real(w) = x + y; Imag(w) = x + y; } // at this point x,y are finite, non-nan. else { // If z is small, then cacos(z) = �/2 - z + z^3/6 + ... = �/2 - z within half an ulp if ((x < 0x1p-13f) && (y < 0x1p-13f)) { Real(w) = pi2f - x; Imag(w) = -y; } // If z is big, then cacos(z) = -i * (log2 + log(z)) + terms smaller than half an ulp else if ((x > 0x1p13f) || (y > 0x1p13f)) { w = clogf(x + I*y) + ln2f; const float tmp = __real__ w; __real__ w = __imag__ w; __imag__ w = -tmp; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Real(w) = 2*atan2( Re(sqrt(1-z)), Re(sqrt(1+z)) ) * Imag(w) = asinh( Im( sqrt(1-z)*sqrt(1+conj(z)) ) ) * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. Analysis of roundoff is a bit harder, but goes though just fine. */ else { ySquared = (y < 0x1p-52f ? 0.0f : y*y); // Avoid underflows. Faster via mask? // Compute sqrt1Plusz = sqrt(1+x + iy) u = 1.0f + x; Real(sqrt1Plusz) = __builtin_sqrtf(0.5f*(__builtin_sqrtf(u*u + ySquared) + u)); Imag(sqrt1Plusz) = 0.5f * (y / Real(sqrt1Plusz)); // Compute sqrt1Minusz = sqrt(1-x - iy) u = 1.0f - x; if (u == 0.0f) { Real(sqrt1Minusz) = sqrt1_2f * __builtin_sqrtf(y); Imag(sqrt1Minusz) = -Real(sqrt1Minusz); } else { Real(sqrt1Minusz) = __builtin_sqrtf(0.5f*(__builtin_sqrtf(u*u + ySquared) + __builtin_fabsf(u))); tmp = 0.5f * (y / Real(sqrt1Minusz)); if (u < 0.0f) { Imag(sqrt1Minusz) = -Real(sqrt1Minusz); Real(sqrt1Minusz) = tmp; } else { Imag(sqrt1Minusz) = -tmp; } } // Magic formulas from Kahan. Real(w) = 2.0f * atan2f(Real(sqrt1Minusz),Real(sqrt1Plusz)); Imag(w) = asinhf( Real(sqrt1Plusz)*Imag(sqrt1Minusz) - Imag(sqrt1Plusz)*Real(sqrt1Minusz) ); } } // Patch up signs to handle z in quadrants II, III & IV, using symmetry. Imag(w) = __builtin_copysignf(Imag(w), -Imag(z)); if (Real(z) < 0.0f) Real(w) = 2.0f * pi2f - Real(w); // No undue cancellation is possible here - Real(w) < �/2. return w;}
long double complex cacosl ( long double complex z ){ static const long double INFl = __builtin_infl(); static const long double ln2l = 0x1.62e42fefa39ef358p-1L; static const long double pi2l = 0x1.921fb54442d1846ap0L; static const long double sqrt1_2l = 0x1.6a09e667f3bcc908p-1L; long double complex w; long double x = __builtin_fabsl(Real(z)); long double y = __builtin_fabsl(Imag(z)); long double u, ySquared, tmp; long double complex sqrt1Plusz, sqrt1Minusz; // If |z| == inf, then cacos(z) = carg(z) - inf i if ((x == INFl) || (y == INFl)) { Imag(w) = -INFl; Real(w) = atan2l(y,x); } // If z = NaN, cacos(z) = NaN, with the special case that cacos(0 + iNaN) = �/2 + iNaN. else if ((x != x) || (y != y)) { if (x == 0.0l) Real(w) = pi2l; else Real(w) = x + y; Imag(w) = x + y; } // at this point x,y are finite, non-nan. else { // If z is small, then cacos(z) = �/2 - z + z^3/6 + ... = �/2 - z within half an ulp if ((x < 0x1p-32l) && (y < 0x1p-32l)) { Real(w) = pi2l - x; Imag(w) = -y; } // If z is big, then cacos(z) = -i * (log2 + log(z)) + terms smaller than half an ulp else if ((x > 0x1p32l) || (y > 0x1p32l)) { w = clogl(x + I*y) + ln2l; const long double tmp = __real__ w; __real__ w = __imag__ w; __imag__ w = -tmp; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Real(w) = 2*atan2( Re(sqrt(1-z)), Re(sqrt(1+z)) ) * Imag(w) = asinh( Im( sqrt(1-z)*sqrt(1+conj(z)) ) ) * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. Analysis of roundoff is a bit harder, but goes though just fine. */ else { ySquared = (y < 0x1p-128l ? 0.0l : y*y); // Avoid underflows. Faster via mask? // Compute sqrt1Plusz = sqrt(1+x + iy) u = 1.0l + x; Real(sqrt1Plusz) = __builtin_sqrtl(0.5l*(__builtin_sqrtl(u*u + ySquared) + u)); Imag(sqrt1Plusz) = 0.5l * (y / Real(sqrt1Plusz)); // Compute sqrt1Minusz = sqrt(1-x - iy) u = 1.0l - x; if (u == 0.0l) { Real(sqrt1Minusz) = sqrt1_2l * __builtin_sqrt(y); Imag(sqrt1Minusz) = -Real(sqrt1Minusz); } else { Real(sqrt1Minusz) = __builtin_sqrtl(0.5l*(__builtin_sqrtl(u*u + ySquared) + __builtin_fabsl(u))); tmp = 0.5l * (y / Real(sqrt1Minusz)); if (u < 0.0l) { Imag(sqrt1Minusz) = -Real(sqrt1Minusz); Real(sqrt1Minusz) = tmp; } else { Imag(sqrt1Minusz) = -tmp; } } // Magic formulas from Kahan. Real(w) = 2.0l * atan2l(Real(sqrt1Minusz), Real(sqrt1Plusz)); Imag(w) = asinhl( Real(sqrt1Plusz)*Imag(sqrt1Minusz) - Imag(sqrt1Plusz)*Real(sqrt1Minusz) ); } } // Patch up signs to handle z in quadrants II, III & IV, using symmetry. Imag(w) = __builtin_copysignl(Imag(w), -Imag(z)); if (Real(z) < 0.0l) Real(w) = 2.0l * pi2l - Real(w); // No undue cancellation is possible here - Real(w) < �/2. return w;}
/**************************************************************************** double complex cacosh(double complex z) returns the complex inverse hyperbolic`cosine of its argument. The algorithm is from Kahan's paper and is based on the formulae: real(cacosh(z)) = arcsinh(real(csqrt(cconj(z)-1.0)*csqrt(z+1.0))) imag(cacosh(z)) = 2.0*atan(imag(csqrt(z-1.0)/real(csqrt(z+1.0)))) Calls: arcsinh, csqrt, atan, feholdexcept, feclearexcept, feupdateenv.****************************************************************************/
double complex cacosh ( double complex z ){ static const double INF = __builtin_inf(); static const double ln2 = 0x1.62e42fefa39ef358p-1; static const double sqrt1_2 = 0x1.6a09e667f3bcc908p-1; static const double pi2 = 0x1.921fb54442d1846ap0; double complex w; double x = __builtin_fabs(Real(z)); double y = __builtin_fabs(Imag(z)); double u, ySquared, tmp; double complex sqrtzPlus1, sqrtzMinus1; // If |z| == inf, then cacosh(z) = inf + carg(z) * I if ((x == INF) || (y == INF)) { Imag(w) = atan2(y,x); Real(w) = INF; } // If z = NaN, cacosh(z) = NaN. else if ((x != x) || (y != y)) { Real(w) = x + y; Imag(w) = x + y; } // at this point x,y are finite, non-nan. else { // If z is small, then cacosh(z) = I*(�/2 - z + z^3/6 + ...) = I*(�/2 - z) within half an ulp if ((x < 0x1p-27) && (y < 0x1p-27)) { Real(w) = y; Imag(w) = pi2 - x; } // If z is big, then cacosh(z) = (log2 + log(z)) + terms smaller than half an ulp else if ((x > 0x1p27) || (y > 0x1p27)) { w = clog(x + I*y) + ln2; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Real(w) = asinh(real(csqrt(cconj(z)-1.0)*csqrt(z+1.0))) * Imag(w) = 2.0*atan2(imag(csqrt(z-1.0))/real(csqrt(z+1.0))) * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. Analysis of roundoff is a bit harder, but goes though just fine. */ else { ySquared = (y < 0x1p-106 ? 0.0 : y*y); // Avoid underflows. Faster via mask? // Compute sqrt1Plusz = sqrt(x+1 + iy) u = x + 1.0; Real(sqrtzPlus1) = __builtin_sqrt(0.5*(__builtin_sqrt(u*u + ySquared) + u)); Imag(sqrtzPlus1) = 0.5 * (y / Real(sqrtzPlus1)); // Compute sqrt1Minusz = sqrt(x-1 + iy) u = x - 1.0; if (u == 0.0) { Real(sqrtzMinus1) = sqrt1_2 * __builtin_sqrt(y); // Avoid spurious underflows in this case Imag(sqrtzMinus1) = Real(sqrtzMinus1); // by using the simpler formula. } else { // No underflow or overflow is possible. Real(sqrtzMinus1) = __builtin_sqrt(0.5*(__builtin_sqrt(u*u + ySquared) + __builtin_fabs(u))); tmp = 0.5 * (y / Real(sqrtzMinus1)); if (u < 0.0) { Imag(sqrtzMinus1) = Real(sqrtzMinus1); Real(sqrtzMinus1) = tmp; } else { Imag(sqrtzMinus1) = tmp; } } // Magic formulas from Kahan. Real(w) = asinh( Real(sqrtzPlus1)*Real(sqrtzMinus1) + Imag(sqrtzPlus1)*Imag(sqrtzMinus1) ); Imag(w) = 2.0*atan2( Imag(sqrtzMinus1) , Real(sqrtzPlus1) ); } } // Patch up signs to handle z in quadrants II, III & IV, using symmetry. if (Real(z) < 0.0) Imag(w) = 2.0 * pi2 - Imag(w); // No undue cancellation is possible here - Real(w) < �/2. Imag(w) = __builtin_copysign(Imag(w), Imag(z)); return w;}
float complex cacoshf ( float complex z ){ static const float INFf = __builtin_inff(); static const float ln2f = 0x1.62e42fefa39ef358p-1f; static const float sqrt1_2f = 0x1.6a09e667f3bcc908p-1f; static const float pi2f = 0x1.921fb54442d1846ap0f; float complex w; float x = __builtin_fabsf(Real(z)); float y = __builtin_fabsf(Imag(z)); float u, ySquared, tmp; float complex sqrtzPlus1, sqrtzMinus1; // If |z| == inf, then cacosh(z) = inf + carg(z) * I if ((x == INFf) || (y == INFf)) { Imag(w) = atan2f(y,x); Real(w) = INFf; } // If z = NaN, cacosh(z) = NaN. else if ((x != x) || (y != y)) { Real(w) = x + y; Imag(w) = x + y; } // at this point x,y are finite, non-nan. else { // If z is small, then cacosh(z) = I*(�/2 - z + z^3/6 + ...) = I*(�/2 - z) within half an ulp if ((x < 0x1p-13f) && (y < 0x1p-13f)) { Real(w) = y; Imag(w) = pi2f - x; } // If z is big, then cacosh(z) = (log2 + log(z)) + terms smaller than half an ulp else if ((x > 0x1p13f) || (y > 0x1p13f)) { w = clogf(x + I*y) + ln2f; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Real(w) = asinh(real(csqrt(cconj(z)-1.0)*csqrt(z+1.0))) * Imag(w) = 2.0*atan2(imag(csqrt(z-1.0))/real(csqrt(z+1.0))) * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. Analysis of roundoff is a bit harder, but goes though just fine. */ else { ySquared = (y < 0x1p-52f ? 0.0f : y*y); // Avoid underflows. Faster via mask? // Compute sqrt1Plusz = sqrt(x+1 + iy) u = x + 1.0f; Real(sqrtzPlus1) = __builtin_sqrtf(0.5f*(__builtin_sqrtf(u*u + ySquared) + u)); Imag(sqrtzPlus1) = 0.5f * (y / Real(sqrtzPlus1)); // Compute sqrt1Minusz = sqrt(x-1 + iy) u = x - 1.0f; if (u == 0.0f) { Real(sqrtzMinus1) = sqrt1_2f * __builtin_sqrtf(y); // Avoid spurious underflows in this case Imag(sqrtzMinus1) = Real(sqrtzMinus1); // by using the simpler formula. } else { // No underflow or overflow is possible. Real(sqrtzMinus1) = __builtin_sqrtf(0.5f*(__builtin_sqrtf(u*u + ySquared) + __builtin_fabsf(u))); tmp = 0.5f * (y / Real(sqrtzMinus1)); if (u < 0.0f) { Imag(sqrtzMinus1) = Real(sqrtzMinus1); Real(sqrtzMinus1) = tmp; } else { Imag(sqrtzMinus1) = tmp; } } // Magic formulas from Kahan. Real(w) = asinhf( Real(sqrtzPlus1)*Real(sqrtzMinus1) + Imag(sqrtzPlus1)*Imag(sqrtzMinus1) ); Imag(w) = 2.0f*atan2f( Imag(sqrtzMinus1) , Real(sqrtzPlus1) ); } } // Patch up signs to handle z in quadrants II, III & IV, using symmetry. if (Real(z) < 0.0f) Imag(w) = 2.0f * pi2f - Imag(w); // No undue cancellation is possible here - Real(w) < �/2. Imag(w) = __builtin_copysignf(Imag(w), Imag(z)); return w;}
long double complex cacoshl ( long double complex z ){ static const long double INFl = __builtin_infl(); static const long double ln2l = 0x1.62e42fefa39ef358p-1L; static const long double sqrt1_2l = 0x1.6a09e667f3bcc908p-1L; static const long double pi2l = 0x1.921fb54442d1846ap0L; long double complex w; long double x = __builtin_fabsl(Real(z)); long double y = __builtin_fabsl(Imag(z)); long double u, ySquared, tmp; long double complex sqrtzPlus1, sqrtzMinus1; // If |z| == inf, then cacosh(z) = inf + carg(z) * I if ((x == INFl) || (y == INFl)) { Imag(w) = atan2l(y,x); Real(w) = INFl; } // If z = NaN, cacosh(z) = NaN. else if ((x != x) || (y != y)) { Real(w) = x + y; Imag(w) = x + y; } // at this point x,y are finite, non-nan. else { // If z is small, then cacosh(z) = I*(�/2 - z + z^3/6 + ...) = I*(�/2 - z) within half an ulp if ((x < 0x1p-32l) && (y < 0x1p-32l)) { Real(w) = y; Imag(w) = pi2l - x; } // If z is big, then cacosh(z) = (log2 + log(z)) + terms smaller than half an ulp else if ((x > 0x1p32l) || (y > 0x1p32l)) { w = clogl(x + I*y) + ln2l; } /* Otherwise, use the expressions from Kahan's "Much ado about nothing's sign bit..." * * Real(w) = asinh(real(csqrt(cconj(z)-1.0)*csqrt(z+1.0))) * Imag(w) = 2.0*atan2(imag(csqrt(z-1.0))/real(csqrt(z+1.0))) * * Derivation of these formulats boggles the mind, but they are easily verified with the * Monodromy theorem. Analysis of roundoff is a bit harder, but goes though just fine. */ else { ySquared = (y < 0x1p-128L ? 0.0L : y*y); // Avoid underflows. Faster via mask? // Compute sqrt1Plusz = sqrt(x+1 + iy) u = x + 1.0l; Real(sqrtzPlus1) = __builtin_sqrtl(0.5l*(__builtin_sqrtl(u*u + ySquared) + u)); Imag(sqrtzPlus1) = 0.5l * (y / Real(sqrtzPlus1)); // Compute sqrt1Minusz = sqrt(x-1 + iy) u = x - 1.0l; if (u == 0.0l) { Real(sqrtzMinus1) = sqrt1_2l * __builtin_sqrtl(y); // Avoid spurious underflows in this case Imag(sqrtzMinus1) = Real(sqrtzMinus1); // by using the simpler formula. } else { // No underflow or overflow is possible. Real(sqrtzMinus1) = __builtin_sqrtl(0.5l*(__builtin_sqrtl(u*u + ySquared) + __builtin_fabsl(u))); tmp = 0.5l * (y / Real(sqrtzMinus1)); if (u < 0.0l) { Imag(sqrtzMinus1) = Real(sqrtzMinus1); Real(sqrtzMinus1) = tmp; } else { Imag(sqrtzMinus1) = tmp; } } // Magic formulas from Kahan. Real(w) = asinhl( Real(sqrtzPlus1)*Real(sqrtzMinus1) + Imag(sqrtzPlus1)*Imag(sqrtzMinus1) ); Imag(w) = 2.0l*atan2l( Imag(sqrtzMinus1) , Real(sqrtzPlus1) ); } } // Patch up signs to handle z in quadrants II, III & IV, using symmetry. if (Real(z) < 0.0l) Imag(w) = 2.0l * pi2l - Imag(w); // No undue cancellation is possible here - Real(w) < �/2. Imag(w) = __builtin_copysignl(Imag(w), Imag(z)); return w;}
/**************************************************************************** double complex catan(double complex z) returns the complex inverse tangent of its argument. The algorithm is from Kahan's paper and is based on the formula: catan(z) = i*(clog(1.0-i*z) - clog(1+i*z))/2.0. CONSTANTS: FPKTHETA = sqrt(nextafterd(+INF,0.0))/4.0 FPKRHO = 1.0/FPKTHETA FPKPI2 = pi/2.0 FPKPI = pi Calls: copysign, fabs, xdivc, sqrt, log, atan, log1p, and carg.****************************************************************************/
double complex catan ( double complex z ){ double complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = catanh(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
float complex catanf ( float complex z ){ float complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = catanhf(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
long double complex catanl ( long double complex z ){ long double complex iz, iw, w; Real(iz) = -Imag(z); Imag(iz) = Real(z); iw = catanhl(iz); Real(w) = Imag(iw); Imag(w) = -Real(iw); return w;}
/**************************************************************************** double complex catanh(double complex z) returns the complex inverse hyperbolic tangent of its argument. The algorithm is from Kahan's paper and is based on the formula: catanh(z) = (clog(1.0 + z) - clog(1 - z))/2.0. CONSTANTS: FPKTHETA = sqrt(nextafterd(+INF,0.0))/4.0 FPKRHO = 1.0/FPKTHETA FPKPI2 = pi/2.0 FPKPI = pi Calls: copysign, fabs, xdivc, sqrt, log, atan, log1p, and carg.****************************************************************************/
double complex catanh( double complex z ) { double complex ctemp, w; double t1, t2, xi, eta, beta; beta = __builtin_copysign(1.0,Real(z)); /* copes with unsigned zero */ Imag(z) = -beta*Imag(z); /* transform real & imag components */ Real(z) = beta*Real(z); if ((Real(z) > FPKTHETA) || (__builtin_fabs(Imag(z)) > FPKTHETA)) { eta = __builtin_copysign(M_PI_2,Imag(z)); /* avoid overflow */ ctemp = xdivc(1.0,z); xi = Real(ctemp); } else if (Real(z) == 1.0) { t1 = __builtin_fabs(Imag(z)) + FPKRHO; xi = log(__builtin_sqrt(__builtin_sqrt(4.0 + t1*t1))/__builtin_sqrt(__builtin_fabs(Imag(z)))); eta = 0.5*__builtin_copysign(M_PI-atan(2.0/(__builtin_fabs(Imag(z))+FPKRHO)),Imag(z)); } else { /* usual case */ t2 = __builtin_fabs(Imag(z)) + FPKRHO; t1 = 1.0 - Real(z); t2 = t2*t2; xi = 0.25*log1p(4.0*Real(z)/(t1*t1 + t2)); Real(ctemp) = (1.0 - Real(z))*(1.0 + Real(z)) - t2; Imag(ctemp) = Imag(z) + Imag(z); eta = 0.5*carg(ctemp); } Real(w) = beta*xi; /* fix up signs of result */ Imag(w) = -beta*eta; return w;}
float complex catanhf( float complex z ) { float complex ctemp, w; float t1, t2, xi, eta, beta; beta = __builtin_copysignf(1.0f,Real(z)); /* copes with unsigned zero */ Imag(z) = -beta*Imag(z); /* transform real & imag components */ Real(z) = beta*Real(z); if ((Real(z) > FPKTHETAf) || (__builtin_fabsf(Imag(z)) > FPKTHETAf)) { eta = __builtin_copysignf((float) M_PI_2,Imag(z)); /* avoid overflow */ ctemp = xdivcf(1.0f,z); xi = Real(ctemp); } else if (Real(z) == 1.0f) { t1 = __builtin_fabsf(Imag(z)) + FPKRHOf; xi = logf(__builtin_sqrtf(__builtin_sqrtf(4.0f + t1*t1))/__builtin_sqrtf(__builtin_fabsf(Imag(z)))); eta = 0.5f*__builtin_copysignf((float)( M_PI-atan(2.0f/(__builtin_fabsf(Imag(z))+FPKRHOf))),Imag(z)); } else { /* usual case */ t2 = __builtin_fabsf(Imag(z)) + FPKRHOf; t1 = 1.0f - Real(z); t2 = t2*t2; xi = 0.25f*log1pf(4.0f*Real(z)/(t1*t1 + t2)); Real(ctemp) = (1.0f - Real(z))*(1.0f + Real(z)) - t2; Imag(ctemp) = Imag(z) + Imag(z); eta = 0.5f*cargf(ctemp); } Real(w) = beta*xi; /* fix up signs of result */ Imag(w) = -beta*eta; return w;}
long double complex catanhl( long double complex z ) { long double complex ctemp, w; long double t1, t2, xi, eta, beta; beta = __builtin_copysignl(1.0L,Real(z)); /* copes with unsigned zero */ Imag(z) = -beta*Imag(z); /* transform real & imag components */ Real(z) = beta*Real(z); if ((Real(z) > FPKTHETA) || (__builtin_fabsl(Imag(z)) > FPKTHETA)) { eta = __builtin_copysignl(M_PI_2,Imag(z)); /* avoid overflow */ ctemp = xdivcl(1.0L,z); xi = Real(ctemp); } else if (Real(z) == 1.0L) { t1 = __builtin_fabsl(Imag(z)) + FPKRHO; xi = logl(__builtin_sqrtl(__builtin_sqrtl(4.0L + t1*t1))/__builtin_sqrtl(__builtin_fabsl(Imag(z)))); eta = 0.5L*__builtin_copysignl(M_PI-atanl(2.0L/(__builtin_fabsl(Imag(z))+FPKRHO)),Imag(z)); } else { /* usual case */ t2 = __builtin_fabsl(Imag(z)) + FPKRHO; t1 = 1.0L - Real(z); t2 = t2*t2; xi = 0.25L*log1pl(4.0L*Real(z)/(t1*t1 + t2)); Real(ctemp) = (1.0L - Real(z))*(1.0L + Real(z)) - t2; Imag(ctemp) = Imag(z) + Imag(z); eta = 0.5L*cargl(ctemp); } Real(w) = beta*xi; /* fix up signs of result */ Imag(w) = -beta*eta; return w;}
/* conj(), creal(), and cimag() are gcc built ins. */double creal( double complex z ){ return __builtin_creal(z);}
float crealf( float complex z ){ return __builtin_crealf(z);}
long double creall( long double complex z ){ return __builtin_creall(z);}
double cimag( double complex z ){ return __builtin_cimag(z);}
float cimagf( float complex z ){ return __builtin_cimagf(z);}
long double cimagl( long double complex z ){ return __builtin_cimagl(z);}
double complex conj( double complex z ){ return __builtin_conj(z);}
float complex conjf( float complex z ){ return __builtin_conjf(z);}
long double complex conjl( long double complex z ){ return __builtin_conjl(z);}
double complex cproj( double complex z ){ static const double inf = __builtin_inf(); double u = __builtin_fabs(Real(z)); double v = __builtin_fabs(Imag(z)); if (EXPECT_FALSE((u == inf) || (v == inf))) { __real__ z = inf; __imag__ z = __builtin_copysign(0.0, __imag__ z); } return z;}
float complex cprojf( float complex z ){ static const float inff = __builtin_inff(); float u = __builtin_fabsf(Real(z)); float v = __builtin_fabsf(Imag(z)); if (EXPECT_FALSE((u == inff) || (v == inff))) { __real__ z = inff; __imag__ z = __builtin_copysignf(0.0f, __imag__ z); } return z;}
long double complex cprojl( long double complex z ){ static const long double infl = __builtin_infl(); long double u = __builtin_fabsl(Real(z)); long double v = __builtin_fabsl(Imag(z)); if (EXPECT_FALSE((u == infl) || (v == infl))) { __real__ z = infl; __imag__ z = __builtin_copysignl(0.0L, __imag__ z); } return z;}