Funkcja wykorzystuje atan2

Kod źródłowy

edytuj

W bibliotece glibc (GNU C Library)[1] jest zdefiniowana w pliku carg.c[2]:

// glibc/math/carg.c
/* Compute argument of complex double value.
   Copyright (C) 1997-2015 Free Software Foundation, Inc.
   This file is part of the GNU C Library.
   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.

   The GNU C Library is free software; you can redistribute it and/or
   modify it under the terms of the GNU Lesser General Public
   License as published by the Free Software Foundation; either
   version 2.1 of the License, or (at your option) any later version.

   The GNU C Library is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with the GNU C Library; if not, see
   <http://www.gnu.org/licenses/>.  */

#include <complex.h>
#include <math.h>

double
__carg (__complex__ double x)
{
  return __atan2 (__imag__ x, __real__ x);
}
weak_alias (__carg, carg)
#ifdef NO_LONG_DOUBLE
strong_alias (__carg, __cargl)
weak_alias (__carg, cargl)
#endif
// http://www.opensource.apple.com/source/Libm/Libm-92/complex.c

/*
 * 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):

****************************************************************************/

/****************************************************************************
   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 )
{
   double a,b,argr;
   int   clre,clim;
   
   a = Real(z);
   b = Imag(z);
   clre = fpclassify(a);
   clim = fpclassify(b);
   
   if ((clre == FP_ZERO) && (clim == FP_ZERO)) {  /* zero real and imag parts */
      a = copysign(1.0, a);
   }
   
   if ((clre == FP_INFINITE) && (clim == FP_INFINITE)) {   /* both parts INF */
      a = copysign(1.0, a);
      b = copysign(1.0, b);
   }
         
   if (fabs(b) > fabs(a))               /* |imag| > |real| */
      argr = copysign(M_PI_2, b) - atan(a/b);

   else {
      if (a < 0.0)                      /* |real| >= |imag| */
         argr = copysign(M_PI, b) + atan(b/a);
      else
         argr = atan(b/a);
   }
   
   return argr;
}

Przykład

edytuj
// gcc c.c -lm -Wall
// a.out

// the return value is the range of [-pi, pi].
// from http://en.cppreference.com/w/c/numeric/complex/carg


#include <stdio.h>
#include <math.h> // M_PI; needs -lm also 
#include <complex.h>

static double TwoPi=2.0*M_PI;

double c_arg(complex double z)
{
 double arg;
 arg = carg(z);
 if (arg<0.0) arg+= TwoPi ; 
 return arg; 
}

double c_turn(complex double z)
{
 double arg;
 arg = c_arg(z);
 return arg/TwoPi; 
}



int info(double complex z)
{
 printf("phase angle of %.1f%+.4fi is %.16f radians  = %0.16f turns\n", creal(z), cimag(z), c_arg(z), c_turn(z));
 return 0;

}


 
int main(void) 
{
    

    info( 1.0+0.0*I);
    //
    info( 1.0+1.0*I);
    //
    info( 0.0+1.0*I);
    //
    info(-1.0+1.0*I);
    //
    info(-1.0+0.0*I);
    info(-1.0-0.0*I);// or CMPLX(-1, -0.0) (the other side of the cut) 
    //
    info( 0.0-1.0*I);
    info( 1.0-1.0*I);
    //
    info( 1.0-0.0001*I);
    info( 1.0-0.0*I);
    
    //
    info( 0.0+0.0*I);


    
 return 0;
}

Wynik :

phase angle of 1.0+0.0000i is 0.0000000000000000 radians  = 0.0000000000000000 turns
phase angle of 1.0+1.0000i is 0.7853981633974483 radians  = 0.1250000000000000 turns
phase angle of 0.0+1.0000i is 1.5707963267948966 radians  = 0.2500000000000000 turns
phase angle of -1.0+1.0000i is 2.3561944901923448 radians  = 0.3750000000000000 turns
phase angle of -1.0+0.0000i is 3.1415926535897931 radians  = 0.5000000000000000 turns
phase angle of -1.0-0.0000i is 3.1415926535897931 radians  = 0.5000000000000000 turns
phase angle of 0.0-1.0000i is 4.7123889803846897 radians  = 0.7500000000000000 turns
phase angle of 1.0-1.0000i is 5.4977871437821380 radians  = 0.8750000000000000 turns
phase angle of 1.0-0.0001i is 6.2830853071799195 radians  = 0.9999840845057438 turns
phase angle of 1.0-0.0000i is -0.0000000000000000 radians  = -0.0000000000000000 turns
phase angle of 0.0+0.0000i is 0.0000000000000000 radians  = 0.0000000000000000 turns


Drugi przykład : obliczamy iMax punktów zespolonych położonych na kole o promieniu r. Nastęþnie sprawdzamy carg

/*



gcc c.c -Wall

./a.out
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>// floor 


double TwoPi=2.0*M_PI;

double GiveTurn(double complex z)
{
  double argument;
  
  argument = carg(z); //   argument in radians from -pi to pi
  
  if ( argument < -M_PI) printf("for z = %f; %f carg argument underflow = %f \n", creal(z), cimag(z), argument);
  if ( argument > M_PI) printf("for z = %f; %f carg argument overflow = %f \n", creal(z), cimag(z), argument);
  
  while (argument<0.0) 
     argument = TwoPi + argument ; //   argument in radians from 0 to 2*pi
     
     
  
  return (argument/TwoPi) ; // argument in turns from 0.0 to 1.0
}



int main(){


int iMax = 30;
int i;
double dt = 1.0/iMax;
double t; // angle in turns
double a; // angle in radians
double tt ; // angle in turns 
double r = 18.0;
double x;
double y;

double arg; 



for (i=0; i<=iMax; i++){

  t = i*dt;
  a = 2.0*t*M_PI; // angle in turns from 0 to 1.0
  // from polar to Rectangular form of complex number z = x+y*I
  x = r* cos(a);
  y = r * sin(a);
  //
  arg = carg(x+y*I);
  
  tt = GiveTurn(x+y*I);
  printf ("i = %d  t = %f tt = %f a = %f  carg = %f z = %f ; %f\n", i, t, tt, a, arg,  x,y );
}



return 0;

}

Wynik:

./a.out
i = 0  t = 0.000000 tt = 0.000000 a = 0.000000  carg = 0.000000 z = 18.000000 ; 0.000000
i = 1  t = 0.033333 tt = 0.033333 a = 0.209440  carg = 0.209440 z = 17.606657 ; 3.742410
i = 2  t = 0.066667 tt = 0.066667 a = 0.418879  carg = 0.418879 z = 16.443818 ; 7.321260
i = 3  t = 0.100000 tt = 0.100000 a = 0.628319  carg = 0.628319 z = 14.562306 ; 10.580135
i = 4  t = 0.133333 tt = 0.133333 a = 0.837758  carg = 0.837758 z = 12.044351 ; 13.376607
i = 5  t = 0.166667 tt = 0.166667 a = 1.047198  carg = 1.047198 z = 9.000000 ; 15.588457
i = 6  t = 0.200000 tt = 0.200000 a = 1.256637  carg = 1.256637 z = 5.562306 ; 17.119017
i = 7  t = 0.233333 tt = 0.233333 a = 1.466077  carg = 1.466077 z = 1.881512 ; 17.901394
i = 8  t = 0.266667 tt = 0.266667 a = 1.675516  carg = 1.675516 z = -1.881512 ; 17.901394
i = 9  t = 0.300000 tt = 0.300000 a = 1.884956  carg = 1.884956 z = -5.562306 ; 17.119017
i = 10  t = 0.333333 tt = 0.333333 a = 2.094395  carg = 2.094395 z = -9.000000 ; 15.588457
i = 11  t = 0.366667 tt = 0.366667 a = 2.303835  carg = 2.303835 z = -12.044351 ; 13.376607
i = 12  t = 0.400000 tt = 0.400000 a = 2.513274  carg = 2.513274 z = -14.562306 ; 10.580135
i = 13  t = 0.433333 tt = 0.433333 a = 2.722714  carg = 2.722714 z = -16.443818 ; 7.321260
i = 14  t = 0.466667 tt = 0.466667 a = 2.932153  carg = 2.932153 z = -17.606657 ; 3.742410
i = 15  t = 0.500000 tt = 0.500000 a = 3.141593  carg = 3.141593 z = -18.000000 ; 0.000000
i = 16  t = 0.533333 tt = 0.533333 a = 3.351032  carg = -2.932153 z = -17.606657 ; -3.742410
i = 17  t = 0.566667 tt = 0.566667 a = 3.560472  carg = -2.722714 z = -16.443818 ; -7.321260
i = 18  t = 0.600000 tt = 0.600000 a = 3.769911  carg = -2.513274 z = -14.562306 ; -10.580135
i = 19  t = 0.633333 tt = 0.633333 a = 3.979351  carg = -2.303835 z = -12.044351 ; -13.376607
i = 20  t = 0.666667 tt = 0.666667 a = 4.188790  carg = -2.094395 z = -9.000000 ; -15.588457
i = 21  t = 0.700000 tt = 0.700000 a = 4.398230  carg = -1.884956 z = -5.562306 ; -17.119017
i = 22  t = 0.733333 tt = 0.733333 a = 4.607669  carg = -1.675516 z = -1.881512 ; -17.901394
i = 23  t = 0.766667 tt = 0.766667 a = 4.817109  carg = -1.466077 z = 1.881512 ; -17.901394
i = 24  t = 0.800000 tt = 0.800000 a = 5.026548  carg = -1.256637 z = 5.562306 ; -17.119017
i = 25  t = 0.833333 tt = 0.833333 a = 5.235988  carg = -1.047198 z = 9.000000 ; -15.588457
i = 26  t = 0.866667 tt = 0.866667 a = 5.445427  carg = -0.837758 z = 12.044351 ; -13.376607
i = 27  t = 0.900000 tt = 0.900000 a = 5.654867  carg = -0.628319 z = 14.562306 ; -10.580135
i = 28  t = 0.933333 tt = 0.933333 a = 5.864306  carg = -0.418879 z = 16.443818 ; -7.321260
i = 29  t = 0.966667 tt = 0.966667 a = 6.073746  carg = -0.209440 z = 17.606657 ; -3.742410
i = 30  t = 1.000000 tt = 1.000000 a = 6.283185  carg = -0.000000 z = 18.000000 ; -0.000000

Źródła

edytuj
  1. glibc download
  2. Kod pliku carg.c