glibc
2.9

00001 00002 /* 00003 * IBM Accurate Mathematical Library 00004 * written by International Business Machines Corp. 00005 * Copyright (C) 2001 Free Software Foundation 00006 * 00007 * This program is free software; you can redistribute it and/or modify 00008 * it under the terms of the GNU Lesser General Public License as published by 00009 * the Free Software Foundation; either version 2.1 of the License, or 00010 * (at your option) any later version. 00011 * 00012 * This program is distributed in the hope that it will be useful, 00013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00015 * GNU Lesser General Public License for more details. 00016 * 00017 * You should have received a copy of the GNU Lesser General Public License 00018 * along with this program; if not, write to the Free Software 00019 * Foundation, Inc., 59 Temple Place  Suite 330, Boston, MA 021111307, USA. 00020 */ 00021 /******************************************************************/ 00022 /* */ 00023 /* MODULE_NAME:mpatan.c */ 00024 /* */ 00025 /* FUNCTIONS:mpatan */ 00026 /* */ 00027 /* FILES NEEDED: mpa.h endian.h mpatan.h */ 00028 /* mpa.c */ 00029 /* */ 00030 /* MultiPrecision Atan function subroutine, for precision p >= 4.*/ 00031 /* The relative error of the result is bounded by 34.32*r**(1p), */ 00032 /* where r=2**24. */ 00033 /******************************************************************/ 00034 00035 #include "endian.h" 00036 #include "mpa.h" 00037 void __mpsqrt(mp_no *, mp_no *, int); 00038 00039 void __mpatan(mp_no *x, mp_no *y, int p) { 00040 #include "mpatan.h" 00041 00042 int i,m,n; 00043 double dx; 00044 mp_no 00045 mpone = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00046 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00047 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, 00048 mptwo = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00049 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00050 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}, 00051 mptwoim1 = {0,{0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00052 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0, 00053 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}}; 00054 00055 mp_no mps,mpsm,mpt,mpt1,mpt2,mpt3; 00056 00057 /* Choose m and initiate mpone, mptwo & mptwoim1 */ 00058 if (EX>0) m=7; 00059 else if (EX<0) m=0; 00060 else { 00061 __mp_dbl(x,&dx,p); dx=ABS(dx); 00062 for (m=6; m>0; m) 00063 {if (dx>xm[m].d) break;} 00064 } 00065 mpone.e = mptwo.e = mptwoim1.e = 1; 00066 mpone.d[0] = mpone.d[1] = mptwo.d[0] = mptwoim1.d[0] = ONE; 00067 mptwo.d[1] = TWO; 00068 00069 /* Reduce x m times */ 00070 __mul(x,x,&mpsm,p); 00071 if (m==0) __cpy(x,&mps,p); 00072 else { 00073 for (i=0; i<m; i++) { 00074 __add(&mpone,&mpsm,&mpt1,p); 00075 __mpsqrt(&mpt1,&mpt2,p); 00076 __add(&mpt2,&mpt2,&mpt1,p); 00077 __add(&mptwo,&mpsm,&mpt2,p); 00078 __add(&mpt1,&mpt2,&mpt3,p); 00079 __dvd(&mpsm,&mpt3,&mpt1,p); 00080 __cpy(&mpt1,&mpsm,p); 00081 } 00082 __mpsqrt(&mpsm,&mps,p); mps.d[0] = X[0]; 00083 } 00084 00085 /* Evaluate a truncated power series for Atan(s) */ 00086 n=np[p]; mptwoim1.d[1] = twonm1[p].d; 00087 __dvd(&mpsm,&mptwoim1,&mpt,p); 00088 for (i=n1; i>1; i) { 00089 mptwoim1.d[1] = TWO; 00090 __dvd(&mpsm,&mptwoim1,&mpt1,p); 00091 __mul(&mpsm,&mpt,&mpt2,p); 00092 __sub(&mpt1,&mpt2,&mpt,p); 00093 } 00094 __mul(&mps,&mpt,&mpt1,p); 00095 __sub(&mps,&mpt1,&mpt,p); 00096 00097 /* Compute Atan(x) */ 00098 mptwoim1.d[1] = twom[m].d; 00099 __mul(&mptwoim1,&mpt,y,p); 00100 00101 return; 00102 }