1/* 2 * ***************************************************************************** 3 * 4 * SPDX-License-Identifier: BSD-2-Clause 5 * 6 * Copyright (c) 2018-2021 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 36scale=20 37define e(x){ 38 auto b,s,n,r,d,i,p,f,v 39 b=ibase 40 ibase=A 41 if(x<0){ 42 n=1 43 x=-x 44 } 45 s=scale 46 r=6+s+.44*x 47 scale=scale(x)+1 48 while(x>1){ 49 d+=1 50 x/=2 51 scale+=1 52 } 53 scale=r 54 r=x+1 55 p=x 56 f=v=1 57 for(i=2;v;++i){ 58 p*=x 59 f*=i 60 v=p/f 61 r+=v 62 } 63 while(d--)r*=r 64 scale=s 65 ibase=b 66 if(n)return(1/r) 67 return(r/1) 68} 69define l(x){ 70 auto b,s,r,p,a,q,i,v 71 if(x<=0)return((1-A^scale)/1) 72 b=ibase 73 ibase=A 74 s=scale 75 scale+=6 76 p=2 77 while(x>=2){ 78 p*=2 79 x=sqrt(x) 80 } 81 while(x<=.5){ 82 p*=2 83 x=sqrt(x) 84 } 85 r=a=(x-1)/(x+1) 86 q=a*a 87 v=1 88 for(i=3;v;i+=2){ 89 a*=q 90 v=a/i 91 r+=v 92 } 93 r*=p 94 scale=s 95 ibase=b 96 return(r/1) 97} 98define s(x){ 99 auto b,s,r,a,q,i 100 if(x<0)return(-s(-x)) 101 b=ibase 102 ibase=A 103 s=scale 104 scale=1.1*s+2 105 a=a(1) 106 scale=0 107 q=(x/a+2)/4 108 x-=4*q*a 109 if(q%2)x=-x 110 scale=s+2 111 r=a=x 112 q=-x*x 113 for(i=3;a;i+=2){ 114 a*=q/(i*(i-1)) 115 r+=a 116 } 117 scale=s 118 ibase=b 119 return(r/1) 120} 121define c(x){ 122 auto b,s 123 b=ibase 124 ibase=A 125 s=scale 126 scale*=1.2 127 x=s(2*a(1)+x) 128 scale=s 129 ibase=b 130 return(x/1) 131} 132define a(x){ 133 auto b,s,r,n,a,m,t,f,i,u 134 b=ibase 135 ibase=A 136 n=1 137 if(x<0){ 138 n=-1 139 x=-x 140 } 141 if(scale<65){ 142 if(x==1){ 143 r=.7853981633974483096156608458198757210492923498437764552437361480/n 144 ibase=b 145 return(r) 146 } 147 if(x==.2){ 148 r=.1973955598498807583700497651947902934475851037878521015176889402/n 149 ibase=b 150 return(r) 151 } 152 } 153 s=scale 154 if(x>.2){ 155 scale+=5 156 a=a(.2) 157 } 158 scale=s+3 159 while(x>.2){ 160 m+=1 161 x=(x-.2)/(1+.2*x) 162 } 163 r=u=x 164 f=-x*x 165 t=1 166 for(i=3;t;i+=2){ 167 u*=f 168 t=u/i 169 r+=t 170 } 171 scale=s 172 ibase=b 173 return((m*a+r)/n) 174} 175define j(n,x){ 176 auto b,s,o,a,i,r,v,f 177 b=ibase 178 ibase=A 179 s=scale 180 scale=0 181 n/=1 182 if(n<0){ 183 n=-n 184 o=n%2 185 } 186 a=1 187 for(i=2;i<=n;++i)a*=i 188 scale=1.5*s 189 a=(x^n)/2^n/a 190 r=v=1 191 f=-x*x/4 192 scale+=length(a)-scale(a) 193 for(i=1;v;++i){ 194 v=v*f/i/(n+i) 195 r+=v 196 } 197 scale=s 198 ibase=b 199 if(o)a=-a 200 return(a*r/1) 201} 202