Home | History | Annotate | Line # | Download | only in complex 
      1 /* $NetBSD: casinf.c,v 1.1 2007/08/20 16:01:31 drochner Exp $ */
      2 
      3 /*-
      4  * Copyright (c) 2007 The NetBSD Foundation, Inc.
      5  * All rights reserved.
      6  *
      7  * This code is derived from software written by Stephen L. Moshier.
      8  * It is redistributed by the NetBSD Foundation by permission of the author.
      9  *
     10  * Redistribution and use in source and binary forms, with or without
     11  * modification, are permitted provided that the following conditions
     12  * are met:
     13  * 1. Redistributions of source code must retain the above copyright
     14  *    notice, this list of conditions and the following disclaimer.
     15  * 2. Redistributions in binary form must reproduce the above copyright
     16  *    notice, this list of conditions and the following disclaimer in the
     17  *    documentation and/or other materials provided with the distribution.
     18  *
     19  * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
     20  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
     21  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
     22  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
     23  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
     24  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
     25  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
     26  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
     27  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
     28  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
     29  * POSSIBILITY OF SUCH DAMAGE.
     30  */
     31 
     32 #include "../src/namespace.h"
     33 #include <complex.h>
     34 #include <math.h>
     35 
     36 #ifdef __weak_alias
     37 __weak_alias(casinf, _casinf)
     38 #endif
     39 
     40 float complex
     41 casinf(float complex z)
     42 {
     43 	float complex w;
     44 	float complex ca, ct, zz, z2;
     45 	float x, y;
     46 
     47 	x = crealf(z);
     48 	y = cimagf(z);
     49 
     50 #if 0 /* MD: test is incorrect, casin(>1) is defined */
     51 	if (y == 0.0f) {
     52 		if (fabsf(x) > 1.0) {
     53 			w = M_PI_2 + 0.0f * I;
     54 #if 0
     55 			mtherr ("casin", DOMAIN);
     56 #endif
     57 		} else {
     58 			w = asinf(x) + 0.0f * I;
     59 		}
     60 		return w;
     61 	}
     62 #endif
     63 
     64 /* Power series expansion */
     65 /*
     66 b = cabsf(z);
     67 if( b < 0.125 )
     68 {
     69 z2.r = (x - y) * (x + y);
     70 z2.i = 2.0 * x * y;
     71 
     72 cn = 1.0;
     73 n = 1.0;
     74 ca.r = x;
     75 ca.i = y;
     76 sum.r = x;
     77 sum.i = y;
     78 do
     79 	{
     80 	ct.r = z2.r * ca.r  -  z2.i * ca.i;
     81 	ct.i = z2.r * ca.i  +  z2.i * ca.r;
     82 	ca.r = ct.r;
     83 	ca.i = ct.i;
     84 
     85 	cn *= n;
     86 	n += 1.0;
     87 	cn /= n;
     88 	n += 1.0;
     89 	b = cn/n;
     90 
     91 	ct.r *= b;
     92 	ct.i *= b;
     93 	sum.r += ct.r;
     94 	sum.i += ct.i;
     95 	b = fabsf(ct.r) + fabsf(ct.i);
     96 	}
     97 while( b > MACHEP );
     98 w->r = sum.r;
     99 w->i = sum.i;
    100 return;
    101 }
    102 */
    103 
    104 
    105 	ca = x + y * I;
    106 	ct = ca * I;
    107 	/* sqrt( 1 - z*z) */
    108 	/* cmul( &ca, &ca, &zz ) */
    109 	/*x * x  -  y * y */
    110 	zz = (x - y) * (x + y) + (2.0f * x * y) * I;
    111 
    112 	zz = 1.0f - crealf(zz) - cimagf(zz) * I;
    113 	z2 = csqrtf(zz);
    114 
    115 	zz = ct + z2;
    116 	zz = clogf(zz);
    117 	/* multiply by 1/i = -i */
    118 	w = zz * (-1.0f * I);
    119 	return w;
    120 }
    121