1/* 2 * ***************************************************************************** 3 * 4 * SPDX-License-Identifier: BSD-2-Clause 5 * 6 * Copyright (c) 2018-2023 Gavin D. Howard and contributors. 7 * 8 * Redistribution and use in source and binary forms, with or without 9 * modification, are permitted provided that the following conditions are met: 10 * 11 * * Redistributions of source code must retain the above copyright notice, this 12 * list of conditions and the following disclaimer. 13 * 14 * * Redistributions in binary form must reproduce the above copyright notice, 15 * this list of conditions and the following disclaimer in the documentation 16 * and/or other materials provided with the distribution. 17 * 18 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 19 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 20 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 21 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 22 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 23 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 24 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 25 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 26 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 27 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 28 * POSSIBILITY OF SUCH DAMAGE. 29 * 30 * ***************************************************************************** 31 * 32 * The bc math library. 33 * 34 */ 35 36define e(x){ 37 auto b,s,n,r,d,i,p,f,v 38 b=ibase 39 ibase=A 40 if(x<0){ 41 n=1 42 x=-x 43 } 44 s=scale 45 r=6+s+.44*x 46 scale=scale(x)+1 47 while(x>1){ 48 d+=1 49 x/=2 50 scale+=1 51 } 52 scale=r 53 r=x+1 54 p=x 55 f=v=1 56 for(i=2;v;++i){ 57 p*=x 58 f*=i 59 v=p/f 60 r+=v 61 } 62 while(d--)r*=r 63 scale=s 64 ibase=b 65 if(n)return(1/r) 66 return(r/1) 67} 68define l(x){ 69 auto b,s,r,p,a,q,i,v 70 if(x<=0)return((1-A^scale)/1) 71 b=ibase 72 ibase=A 73 s=scale 74 scale+=6 75 p=2 76 while(x>=2){ 77 p*=2 78 x=sqrt(x) 79 } 80 while(x<=.5){ 81 p*=2 82 x=sqrt(x) 83 } 84 r=a=(x-1)/(x+1) 85 q=a*a 86 v=1 87 for(i=3;v;i+=2){ 88 a*=q 89 v=a/i 90 r+=v 91 } 92 r*=p 93 scale=s 94 ibase=b 95 return(r/1) 96} 97define s(x){ 98 auto b,s,r,a,q,i 99 if(x<0)return(-s(-x)) 100 b=ibase 101 ibase=A 102 s=scale 103 scale=1.1*s+2 104 a=a(1) 105 scale=0 106 q=(x/a+2)/4 107 x-=4*q*a 108 if(q%2)x=-x 109 scale=s+2 110 r=a=x 111 q=-x*x 112 for(i=3;a;i+=2){ 113 a*=q/(i*(i-1)) 114 r+=a 115 } 116 scale=s 117 ibase=b 118 return(r/1) 119} 120define c(x){ 121 auto b,s 122 b=ibase 123 ibase=A 124 s=scale 125 scale*=1.2 126 x=s(2*a(1)+x) 127 scale=s 128 ibase=b 129 return(x/1) 130} 131define a(x){ 132 auto b,s,r,n,a,m,t,f,i,u 133 b=ibase 134 ibase=A 135 n=1 136 if(x<0){ 137 n=-1 138 x=-x 139 } 140 if(scale<65){ 141 if(x==1){ 142 r=.7853981633974483096156608458198757210492923498437764552437361480/n 143 ibase=b 144 return(r) 145 } 146 if(x==.2){ 147 r=.1973955598498807583700497651947902934475851037878521015176889402/n 148 ibase=b 149 return(r) 150 } 151 } 152 s=scale 153 if(x>.2){ 154 scale+=5 155 a=a(.2) 156 } 157 scale=s+3 158 while(x>.2){ 159 m+=1 160 x=(x-.2)/(1+.2*x) 161 } 162 r=u=x 163 f=-x*x 164 t=1 165 for(i=3;t;i+=2){ 166 u*=f 167 t=u/i 168 r+=t 169 } 170 scale=s 171 ibase=b 172 return((m*a+r)/n) 173} 174define j(n,x){ 175 auto b,s,o,a,i,r,v,f 176 b=ibase 177 ibase=A 178 s=scale 179 scale=0 180 n/=1 181 if(n<0){ 182 n=-n 183 o=n%2 184 } 185 a=1 186 for(i=2;i<=n;++i)a*=i 187 scale=1.5*s 188 a=(x^n)/2^n/a 189 r=v=1 190 f=-x*x/4 191 scale+=length(a)-scale(a) 192 for(i=1;v;++i){ 193 v=v*f/i/(n+i) 194 r+=v 195 } 196 scale=s 197 ibase=b 198 if(o)a=-a 199 return(a*r/1) 200} 201