Home | History | Annotate | Download | only in bc
      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 
     32 scale = 20
     33 define 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 
     79 define 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 
    129 define 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 
    160 define 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 
    172 define 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 
    213 define j(n,x){
    214 auto 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 
    223 s= -x*x/4
    224 if(n<0){
    225 	n= -n
    226 	x= -x
    227 	}
    228 a=1
    229 c=1
    230 for(i=1;i<=n;i++){
    231 	a=a*x
    232 	c = c*2*i
    233 	}
    234 b=a
    235 d=1
    236 e=1
    237 for(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