← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to deps/openlibm/slatec/xpqnu.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-01-16 12:29:42 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130116122942-x86e42akjq31repw
Tags: 0.0.0+20130107.gitd9656f41-1
* New upstream snashot
* No longer try to rebuild helpdb.jl.
   + debian/rules: remove helpdb.jl from build-arch rule
   + debian/control: move back python-sphinx to Build-Depends-Indep
* debian/copyright: reflect upstream changes
* Add Build-Conflicts on libatlas3-base (makes linalg tests fail)
* debian/rules: replace obsolete USE_DEBIAN makeflag by a list of
  USE_SYSTEM_* flags
* debian/rules: on non-x86 systems, use libm instead of openlibm
* dpkg-buildflags.patch: remove patch, applied upstream
* Refreshed other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
deps/openlibm/slatec/xpqnu.f
 
1
*DECK XPQNU
 
2
      SUBROUTINE XPQNU (NU1, NU2, MU, THETA, ID, PQA, IPQA, IERROR)
 
3
C***BEGIN PROLOGUE  XPQNU
 
4
C***SUBSIDIARY
 
5
C***PURPOSE  To compute the values of Legendre functions for XLEGF.
 
6
C            This subroutine calculates initial values of P or Q using
 
7
C            power series, then performs forward nu-wise recurrence to
 
8
C            obtain P(-MU,NU,X), Q(0,NU,X), or Q(1,NU,X). The nu-wise
 
9
C            recurrence is stable for P for all mu and for Q for mu=0,1.
 
10
C***LIBRARY   SLATEC
 
11
C***CATEGORY  C3A2, C9
 
12
C***TYPE      SINGLE PRECISION (XPQNU-S, DXPQNU-D)
 
13
C***KEYWORDS  LEGENDRE FUNCTIONS
 
14
C***AUTHOR  Smith, John M., (NBS and George Mason University)
 
15
C***ROUTINES CALLED  XADD, XADJ, XPSI
 
16
C***COMMON BLOCKS    XBLK1
 
17
C***REVISION HISTORY  (YYMMDD)
 
18
C   820728  DATE WRITTEN
 
19
C   890126  Revised to meet SLATEC CML recommendations.  (DWL and JMS)
 
20
C   901019  Revisions to prologue.  (DWL and WRB)
 
21
C   901106  Changed all specific intrinsics to generic.  (WRB)
 
22
C           Corrected order of sections in prologue and added TYPE
 
23
C           section.  (WRB)
 
24
C   920127  Revised PURPOSE section of prologue.  (DWL)
 
25
C***END PROLOGUE  XPQNU
 
26
      REAL A,NU,NU1,NU2,PQ,PQA,XPSI,R,THETA,W,X,X1,X2,XS,
 
27
     1 Y,Z
 
28
      REAL DI,DMU,PQ1,PQ2,FACTMU,FLOK
 
29
      DIMENSION PQA(*),IPQA(*)
 
30
      COMMON /XBLK1/ NBITSF
 
31
      SAVE /XBLK1/
 
32
C
 
33
C        J0, IPSIK, AND IPSIX ARE INITIALIZED IN THIS SUBROUTINE.
 
34
C        J0 IS THE NUMBER OF TERMS USED IN SERIES EXPANSION
 
35
C        IN SUBROUTINE XPQNU.
 
36
C        IPSIK, IPSIX ARE VALUES OF K AND X RESPECTIVELY
 
37
C        USED IN THE CALCULATION OF THE XPSI FUNCTION.
 
38
C
 
39
C***FIRST EXECUTABLE STATEMENT  XPQNU
 
40
      IERROR=0
 
41
      J0=NBITSF
 
42
      IPSIK=1+(NBITSF/10)
 
43
      IPSIX=5*IPSIK
 
44
      IPQ=0
 
45
C        FIND NU IN INTERVAL [-.5,.5) IF ID=2  ( CALCULATION OF Q )
 
46
      NU=MOD(NU1,1.)
 
47
      IF(NU.GE..5) NU=NU-1.
 
48
C        FIND NU IN INTERVAL (-1.5,-.5] IF ID=1,3, OR 4  ( CALC. OF P )
 
49
      IF(ID.NE.2.AND.NU.GT.-.5) NU=NU-1.
 
50
C        CALCULATE MU FACTORIAL
 
51
      K=MU
 
52
      DMU=MU
 
53
      IF(MU.LE.0) GO TO 60
 
54
      FACTMU=1.
 
55
      IF=0
 
56
      DO 50 I=1,K
 
57
      FACTMU=FACTMU*I
 
58
   50 CALL XADJ(FACTMU,IF,IERROR)
 
59
      IF (IERROR.NE.0) RETURN
 
60
   60 IF(K.EQ.0) FACTMU=1.
 
61
      IF(K.EQ.0) IF=0
 
62
C
 
63
C        X=COS(THETA)
 
64
C        Y=SIN(THETA/2)**2=(1-X)/2=.5-.5*X
 
65
C        R=TAN(THETA/2)=SQRT((1-X)/(1+X)
 
66
C
 
67
      X=COS(THETA)
 
68
      Y=SIN(THETA/2.)**2
 
69
      R=TAN(THETA/2.)
 
70
C
 
71
C        USE ASCENDING SERIES TO CALCULATE TWO VALUES OF P OR Q
 
72
C        FOR USE AS STARTING VALUES IN RECURRENCE RELATION.
 
73
C
 
74
      PQ2=0.0
 
75
      DO 100 J=1,2
 
76
      IPQ1=0
 
77
      IF(ID.EQ.2) GO TO 80
 
78
C
 
79
C        SERIES FOR P ( ID = 1, 3, OR 4 )
 
80
C        P(-MU,NU,X)=1./FACTORIAL(MU)*SQRT(((1.-X)/(1.+X))**MU)
 
81
C                *SUM(FROM 0 TO J0-1)A(J)*(.5-.5*X)**J
 
82
C
 
83
      IPQ=0
 
84
      PQ=1.
 
85
      A=1.
 
86
      IA=0
 
87
      DO 65 I=2,J0
 
88
      DI=I
 
89
      A=A*Y*(DI-2.-NU)*(DI-1.+NU)/((DI-1.+DMU)*(DI-1.))
 
90
      CALL XADJ(A,IA,IERROR)
 
91
      IF (IERROR.NE.0) RETURN
 
92
      IF(A.EQ.0.) GO TO 66
 
93
      CALL XADD(PQ,IPQ,A,IA,PQ,IPQ,IERROR)
 
94
      IF (IERROR.NE.0) RETURN
 
95
   65 CONTINUE
 
96
   66 CONTINUE
 
97
      IF(MU.LE.0) GO TO 90
 
98
      X2=R
 
99
      X1=PQ
 
100
      K=MU
 
101
      DO 77 I=1,K
 
102
      X1=X1*X2
 
103
   77 CALL XADJ(X1,IPQ,IERROR)
 
104
      IF (IERROR.NE.0) RETURN
 
105
      PQ=X1/FACTMU
 
106
      IPQ=IPQ-IF
 
107
      CALL XADJ(PQ,IPQ,IERROR)
 
108
      IF (IERROR.NE.0) RETURN
 
109
      GO TO 90
 
110
C
 
111
C        Z=-LN(R)=.5*LN((1+X)/(1-X))
 
112
C
 
113
   80 Z=-LOG(R)
 
114
      W=XPSI(NU+1.,IPSIK,IPSIX)
 
115
      XS=1./SIN(THETA)
 
116
C
 
117
C        SERIES SUMMATION FOR Q ( ID = 2 )
 
118
C        Q(0,NU,X)=SUM(FROM 0 TO J0-1)((.5*LN((1+X)/(1-X))
 
119
C    +XPSI(J+1,IPSIK,IPSIX)-XPSI(NU+1,IPSIK,IPSIX)))*A(J)*(.5-.5*X)**J
 
120
C
 
121
C        Q(1,NU,X)=-SQRT(1./(1.-X**2))+SQRT((1-X)/(1+X))
 
122
C             *SUM(FROM 0 T0 J0-1)(-NU*(NU+1)/2*LN((1+X)/(1-X))
 
123
C                 +(J-NU)*(J+NU+1)/(2*(J+1))+NU*(NU+1)*
 
124
C     (XPSI(NU+1,IPSIK,IPSIX)-XPSI(J+1,IPSIK,IPSIX))*A(J)*(.5-.5*X)**J
 
125
C
 
126
C        NOTE, IN THIS LOOP K=J+1
 
127
C
 
128
      PQ=0.
 
129
      IPQ=0
 
130
      IA=0
 
131
      A=1.
 
132
      DO 85 K=1,J0
 
133
      FLOK=K
 
134
      IF(K.EQ.1) GO TO 81
 
135
      A=A*Y*(FLOK-2.-NU)*(FLOK-1.+NU)/((FLOK-1.+DMU)*(FLOK-1.))
 
136
      CALL XADJ(A,IA,IERROR)
 
137
      IF (IERROR.NE.0) RETURN
 
138
   81 CONTINUE
 
139
      IF(MU.GE.1) GO TO 83
 
140
      X1=(XPSI(FLOK,IPSIK,IPSIX)-W+Z)*A
 
141
      IX1=IA
 
142
      CALL XADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
 
143
      IF (IERROR.NE.0) RETURN
 
144
      GO TO 85
 
145
   83 X1=(NU*(NU+1.)*(Z-W+XPSI(FLOK,IPSIK,IPSIX))+(NU-FLOK+1.)
 
146
     1  *(NU+FLOK)/(2.*FLOK))*A
 
147
      IX1=IA
 
148
      CALL XADD(PQ,IPQ,X1,IX1,PQ,IPQ,IERROR)
 
149
      IF (IERROR.NE.0) RETURN
 
150
   85 CONTINUE
 
151
      IF(MU.GE.1) PQ=-R*PQ
 
152
      IXS=0
 
153
      IF(MU.GE.1) CALL XADD(PQ,IPQ,-XS,IXS,PQ,IPQ,IERROR)
 
154
      IF (IERROR.NE.0) RETURN
 
155
      IF(J.EQ.2) MU=-MU
 
156
      IF(J.EQ.2) DMU=-DMU
 
157
   90 IF(J.EQ.1) PQ2=PQ
 
158
      IF(J.EQ.1) IPQ2=IPQ
 
159
      NU=NU+1.
 
160
  100 CONTINUE
 
161
      K=0
 
162
      IF(NU-1.5.LT.NU1) GO TO 120
 
163
      K=K+1
 
164
      PQA(K)=PQ2
 
165
      IPQA(K)=IPQ2
 
166
      IF(NU.GT.NU2+.5) RETURN
 
167
  120 PQ1=PQ
 
168
      IPQ1=IPQ
 
169
      IF(NU.LT.NU1+.5) GO TO 130
 
170
      K=K+1
 
171
      PQA(K)=PQ
 
172
      IPQA(K)=IPQ
 
173
      IF(NU.GT.NU2+.5) RETURN
 
174
C
 
175
C        FORWARD NU-WISE RECURRENCE FOR F(MU,NU,X) FOR FIXED MU
 
176
C        USING
 
177
C        (NU+MU+1)*F(MU,NU,X)=(2.*NU+1)*F(MU,NU,X)-(NU-MU)*F(MU,NU-1,X)
 
178
C        WHERE F(MU,NU,X) MAY BE P(-MU,NU,X) OR IF MU IS REPLACED
 
179
C        BY -MU THEN F(MU,NU,X) MAY BE Q(MU,NU,X).
 
180
C        NOTE, IN THIS LOOP, NU=NU+1
 
181
C
 
182
  130 X1=(2.*NU-1.)/(NU+DMU)*X*PQ1
 
183
      X2=(NU-1.-DMU)/(NU+DMU)*PQ2
 
184
      CALL XADD(X1,IPQ1,-X2,IPQ2,PQ,IPQ,IERROR)
 
185
      IF (IERROR.NE.0) RETURN
 
186
      CALL XADJ(PQ,IPQ,IERROR)
 
187
      IF (IERROR.NE.0) RETURN
 
188
      NU=NU+1.
 
189
      PQ2=PQ1
 
190
      IPQ2=IPQ1
 
191
      GO TO 120
 
192
C
 
193
      END