1/* 2 * ***************************************************************************** 3 * 4 * Copyright (c) 2018-2019 Gavin D. Howard and contributors. 5 * 6 * All rights reserved. 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 second bc math library. 33 * 34 */ 35 36define r(x,p){ 37 auto t,n 38 if(x==0)return x 39 p=abs(p)$ 40 n=(x<0) 41 x=abs(x) 42 t=x@p 43 if(p<scale(x)&&x-t>=5>>p+1)t+=1>>p 44 if(n)t=-t 45 return t 46} 47define ceil(x,p){ 48 auto t,n 49 if(x==0)return x 50 p=abs(p)$ 51 n=(x<0) 52 x=abs(x) 53 t=(x+(9>>p+1))@p 54 if(n)t=-t 55 return t 56} 57define f(n){ 58 auto r 59 n=abs(n)$ 60 for(r=1;n>1;--n)r*=n 61 return r 62} 63define perm(n,k){ 64 auto f,g,s 65 if(k>n)return 0 66 n=abs(n)$ 67 k=abs(k)$ 68 f=f(n) 69 g=f(n-k) 70 s=scale 71 scale=0 72 f/=g 73 scale=s 74 return f 75} 76define comb(n,r){ 77 auto s,f,g,h 78 if(r>n)return 0 79 n=abs(n)$ 80 r=abs(r)$ 81 s=scale 82 scale=0 83 f=f(n) 84 h=f(r) 85 g=f(n-r) 86 f/=h*g 87 scale=s 88 return f 89} 90define log(x,b){ 91 auto p,s 92 s=scale 93 if(scale<K)scale=K 94 if(scale(x)>scale)scale=scale(x) 95 scale*=2 96 p=l(x)/l(b) 97 scale=s 98 return p 99} 100define l2(x){return log(x,2)} 101define l10(x){return log(x,A)} 102define root(x,n){ 103 auto s,m,r,q,p 104 if(n<0)sqrt(n) 105 n=abs(n)$ 106 if(n==0)x/n 107 if(n==1)return x 108 if(n==2)return sqrt(x) 109 s=scale 110 scale=0 111 if(x<0&&n%2==0)sqrt(x) 112 scale=s+2 113 m=(x<0) 114 x=abs(x) 115 p=n-1 116 q=10^ceil((length(x$)/n)$,0) 117 while(r!=q){ 118 r=q 119 q=(p*r+x/r^p)/n 120 } 121 if(m)r=-r 122 scale=s 123 return r@s 124} 125define cbrt(x){return root(x,3)} 126define pi(s){ 127 auto t,v 128 if(s==0)return 3 129 s=abs(s)$ 130 t=scale 131 scale=s+1 132 v=4*a(1) 133 scale=t 134 return v@s 135} 136define t(x){ 137 auto s,c,l 138 l=scale 139 scale+=2 140 s=s(x) 141 c=c(x) 142 scale=l 143 return s/c 144} 145define a2(y,x){ 146 auto a,p 147 if(!x&&!y)y/x 148 if(x<=0){ 149 p=pi(scale+2) 150 if(y<0)p=-p 151 } 152 if(x==0)a=p/2 153 else{ 154 scale+=2 155 a=a(y/x)+p 156 scale-=2 157 } 158 return a@scale 159} 160define sin(x){return s(x)} 161define cos(x){return c(x)} 162define atan(x){return a(x)} 163define tan(x){return t(x)} 164define atan2(y,x){return a2(y,x)} 165define r2d(x){ 166 auto r,i,s 167 s=scale 168 scale+=5 169 i=ibase 170 ibase=A 171 r=x*180/pi(scale) 172 ibase=i 173 scale=s 174 return r@s 175} 176define d2r(x){ 177 auto r,i,s 178 s=scale 179 scale+=5 180 i=ibase 181 ibase=A 182 r=x*pi(scale)/180 183 ibase=i 184 scale=s 185 return r@s 186} 187define void output(x,b){ 188 auto c 189 c=obase 190 obase=b 191 x 192 obase=c 193} 194define void hex(x){output(x,G)} 195define void binary(x){output(x,2)} 196define ubytes(x){ 197 auto p,b,i 198 b=ibase 199 ibase=A 200 x=abs(x)$ 201 i=2^8 202 for(p=1;i-1<x;p*=2){i*=i} 203 ibase=b 204 return p 205} 206define sbytes(x){ 207 auto p,b,n,z 208 z=(x<0) 209 x=abs(x) 210 x=x$ 211 n=ubytes(x) 212 b=ibase 213 ibase=A 214 p=2^(n*8-1) 215 if(x>p||(!z&&x==p))n*=2 216 ibase=b 217 return n 218} 219define void output_byte(x,i){ 220 auto j,p,y,b 221 j=ibase 222 ibase=A 223 s=scale 224 scale=0 225 x=abs(x)$ 226 b=x/(2^(i*8)) 227 b%=2^8 228 y=log(256,obase)$ 229 if(b>1)p=log(b,obase)$+1 230 else p=b 231 for(i=y-p;i>0;--i)print 0 232 if(b)print b 233 scale=s 234 ibase=j 235} 236define void output_uint(x,n){ 237 auto i,b 238 b=ibase 239 ibase=A 240 for(i=n-1;i>=0;--i){ 241 output_byte(x,i) 242 if(i)print" " 243 else print"\n" 244 } 245 ibase=b 246} 247define void hex_uint(x,n){ 248 auto o 249 o=obase 250 obase=G 251 output_uint(x,n) 252 obase=o 253} 254define void binary_uint(x,n){ 255 auto o 256 o=obase 257 obase=2 258 output_uint(x,n) 259 obase=o 260} 261define void uintn(x,n){ 262 if(scale(x)){ 263 print"Error: ",x," is not an integer.\n" 264 return 265 } 266 if(x<0){ 267 print"Error: ",x," is negative.\n" 268 return 269 } 270 if(x>=2^(n*8)){ 271 print"Error: ",x," cannot fit into ",n," unsigned byte(s).\n" 272 return 273 } 274 binary_uint(x,n) 275 hex_uint(x,n) 276} 277define void intn(x,n){ 278 auto t 279 if(scale(x)){ 280 print"Error: ",x," is not an integer.\n" 281 return 282 } 283 t=2^(n*8-1) 284 if(abs(x)>=t&&(x>0||x!=-t)){ 285 print "Error: ",x," cannot fit into ",n," signed byte(s).\n" 286 return 287 } 288 if(x<0)x=2^(n*8)-(-x) 289 binary_uint(x,n) 290 hex_uint(x,n) 291} 292define void uint8(x){uintn(x,1)} 293define void int8(x){intn(x,1)} 294define void uint16(x){uintn(x,2)} 295define void int16(x){intn(x,2)} 296define void uint32(x){uintn(x,4)} 297define void int32(x){intn(x,4)} 298define void uint64(x){uintn(x,8)} 299define void int64(x){intn(x,8)} 300define void uint(x){uintn(x,ubytes(x))} 301define void int(x){intn(x,sbytes(x))} 302