1/* Copyright (c) 1984 AT&T */ 2/* All Rights Reserved */ 3 4 5/* 6 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 7 * Use is subject to license terms. 8 * 9 * CDDL HEADER START 10 * 11 * The contents of this file are subject to the terms of the 12 * Common Development and Distribution License, Version 1.0 only 13 * (the "License"). You may not use this file except in compliance 14 * with the License. 15 * 16 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 17 * or http://www.opensolaris.org/os/licensing. 18 * See the License for the specific language governing permissions 19 * and limitations under the License. 20 * 21 * When distributing Covered Code, include this CDDL HEADER in each 22 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 23 * If applicable, add the following below this CDDL HEADER, with the 24 * fields enclosed by brackets "[]" replaced with your own identifying 25 * information: Portions Copyright [yyyy] [name of copyright owner] 26 * 27 * CDDL HEADER END 28 */ 29 30/* #ident "%Z%%M% %I% %E% SMI" */ 31 32scale = 20 33define e(x){ 34 auto a, b, c, d, e, g, w, y, t, r 35 r = ibase 36 ibase = A 37 38 t = scale 39 scale = t + .434*x + 1 40 41 w = 0 42 if(x<0){ 43 x = -x 44 w = 1 45 } 46 y = 0 47 while(x>2){ 48 x = x/2 49 y = y + 1 50 } 51 52 a=1 53 b=1 54 c=b 55 d=1 56 e=1 57 for(a=1;1==1;a++){ 58 b=b*x 59 c=c*a+b 60 d=d*a 61 g = c/d 62 if(g == e){ 63 g = g/1 64 while(y--){ 65 g = g*g 66 } 67 scale = t 68 if(w==1){ 69 ibase = r 70 return(1/g) 71 } 72 ibase = r 73 return(g/1) 74 } 75 e=g 76 } 77} 78 79define l(x){ 80 auto a, b, c, d, e, f, g, u, s, t, r, z 81 r = ibase 82 ibase = A 83 if(x <=0){ 84 z = 1-10^scale 85 ibase = r 86 return(z) 87 } 88 t = scale 89 90 f=1 91 scale = scale + scale(x) - length(x) + 1 92 s=scale 93 while(x > 2){ 94 s = s + (length(x)-scale(x))/2 + 1 95 if(s>0) scale = s 96 x = sqrt(x) 97 f=f*2 98 } 99 while(x < .5){ 100 s = s + (length(x)-scale(x))/2 + 1 101 if(s>0) scale = s 102 x = sqrt(x) 103 f=f*2 104 } 105 106 scale = t + length(f) - scale(f) + 1 107 u = (x-1)/(x+1) 108 109 scale = scale + 1.1*length(t) - 1.1*scale(t) 110 s = u*u 111 b = 2*f 112 c = b 113 d = 1 114 e = 1 115 for(a=3;1==1;a=a+2){ 116 b=b*s 117 c=c*a+d*b 118 d=d*a 119 g=c/d 120 if(g==e){ 121 scale = t 122 ibase = r 123 return(u*c/d) 124 } 125 e=g 126 } 127} 128 129define s(x){ 130 auto a, b, c, s, t, y, p, n, i, r 131 r = ibase 132 ibase = A 133 t = scale 134 y = x/.7853 135 s = t + length(y) - scale(y) 136 if(s<t) s=t 137 scale = s 138 p = a(1) 139 140 scale = 0 141 if(x>=0) n = (x/(2*p)+1)/2 142 if(x<0) n = (x/(2*p)-1)/2 143 x = x - 4*n*p 144 if(n%2!=0) x = -x 145 146 scale = t + length(1.2*t) - scale(1.2*t) 147 y = -x*x 148 a = x 149 b = 1 150 s = x 151 for(i=3; 1==1; i=i+2){ 152 a = a*y 153 b = b*i*(i-1) 154 c = a/b 155 if(c==0){scale=t; ibase = r; return(s/1)} 156 s = s+c 157 } 158} 159 160define c(x){ 161 auto t, r 162 r = ibase 163 ibase = A 164 t = scale 165 scale = scale+1 166 x = s(x+2*a(1)) 167 scale = t 168 ibase = r 169 return(x/1) 170} 171 172define a(x){ 173 auto a, b, c, d, e, f, g, s, t, r, z 174 if(x==0) {return(0/1)} 175 r = ibase 176 ibase = A 177 if(x==1){ 178 z =.7853981633974483096156608458198757210492923498437764/1 179 ibase = r 180 if(scale<52)return(z) 181 } 182 t = scale 183 f=1 184 while(x > .5){ 185 scale = scale + 1 186 x= -(1-sqrt(1.+x*x))/x 187 f=f*2 188 } 189 while(x < -.5){ 190 scale = scale + 1 191 x = -(1-sqrt(1.+x*x))/x 192 f=f*2 193 } 194 s = -x*x 195 b = f 196 c = f 197 d = 1 198 e = 1 199 for(a=3;1==1;a=a+2){ 200 b=b*s 201 c=c*a+d*b 202 d=d*a 203 g=c/d 204 if(g==e){ 205 scale = t 206 ibase = r 207 return(x*c/d) 208 } 209 e=g 210 } 211} 212 213define j(n,x){ 214auto a,b,c,d,e,g,i,s,k,t, r 215 r = ibase 216 ibase = A 217 218 t = scale 219 k = 1.36*x + 1.16*t - n 220 k = length(k) - scale(k) 221 if(k>0) scale = scale + k 222 223s= -x*x/4 224if(n<0){ 225 n= -n 226 x= -x 227 } 228a=1 229c=1 230for(i=1;i<=n;i++){ 231 a=a*x 232 c = c*2*i 233 } 234b=a 235d=1 236e=1 237for(i=1;1;i++){ 238 a=a*s 239 b=b*i*(n+i) + a 240 c=c*i*(n+i) 241 g=b/c 242 if(g==e){ 243 scale = t 244 ibase = r 245 return(g/1) 246 } 247 e=g 248 } 249} 250 251