contrib/bc/gen/lib2.bc

   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 second bc math library.
  33  *
  34  */
  35
  36 define p(x,y){
  37         auto a
  38         a=y$
  39         if(y==a)return (x^a)@scale
  40         return e(y*l(x))
  41 }
  42 define r(x,p){
  43         auto t,n
  44         if(x==0)return x
  45         p=abs(p)$
  46         n=(x<0)
  47         x=abs(x)
  48         t=x@p
  49         if(p<scale(x)&&x-t>=5>>p+1)t+=1>>p
  50         if(n)t=-t
  51         return t
  52 }
  53 define ceil(x,p){
  54         auto t,n
  55         if(x==0)return x
  56         p=abs(p)$
  57         n=(x<0)
  58         x=abs(x)
  59         t=(x+((x@p<x)>>p))@p
  60         if(n)t=-t
  61         return t
  62 }
  63 define f(n){
  64         auto r
  65         n=abs(n)$
  66         for(r=1;n>1;--n)r*=n
  67         return r
  68 }
  69 define perm(n,k){
  70         auto f,g,s
  71         if(k>n)return 0
  72         n=abs(n)$
  73         k=abs(k)$
  74         f=f(n)
  75         g=f(n-k)
  76         s=scale
  77         scale=0
  78         f/=g
  79         scale=s
  80         return f
  81 }
  82 define comb(n,r){
  83         auto s,f,g,h
  84         if(r>n)return 0
  85         n=abs(n)$
  86         r=abs(r)$
  87         s=scale
  88         scale=0
  89         f=f(n)
  90         h=f(r)
  91         g=f(n-r)
  92         f/=h*g
  93         scale=s
  94         return f
  95 }
  96 define log(x,b){
  97         auto p,s
  98         s=scale
  99         if(scale<K)scale=K
 100         if(scale(x)>scale)scale=scale(x)
 101         scale*=2
 102         p=l(x)/l(b)
 103         scale=s
 104         return p@s
 105 }
 106 define l2(x){return log(x,2)}
 107 define l10(x){return log(x,A)}
 108 define root(x,n){
 109         auto s,m,r,q,p
 110         if(n<0)sqrt(n)
 111         n=n$
 112         if(n==0)x/n
 113         if(x==0||n==1)return x
 114         if(n==2)return sqrt(x)
 115         s=scale
 116         scale=0
 117         if(x<0&&n%2==0)sqrt(x)
 118         scale=s+2
 119         m=(x<0)
 120         x=abs(x)
 121         p=n-1
 122         q=A^ceil((length(x$)/n)$,0)
 123         while(r!=q){
 124                 r=q
 125                 q=(p*r+x/r^p)/n
 126         }
 127         if(m)r=-r
 128         scale=s
 129         return r@s
 130 }
 131 define cbrt(x){return root(x,3)}
 132 define gcd(a,b){
 133         auto g,s
 134         if(!b)return a
 135         s=scale
 136         scale=0
 137         a=abs(a)$
 138         b=abs(b)$
 139         if(a<b){
 140                 g=a
 141                 a=b
 142                 b=g
 143         }
 144         while(b){
 145                 g=a%b
 146                 a=b
 147                 b=g
 148         }
 149         scale=s
 150         return a
 151 }
 152 define lcm(a,b){
 153         auto r,s
 154         if(!a&&!b)return 0
 155         s=scale
 156         scale=0
 157         a=abs(a)$
 158         b=abs(b)$
 159         r=a*b/gcd(a,b)
 160         scale=s
 161         return r
 162 }
 163 define pi(s){
 164         auto t,v
 165         if(s==0)return 3
 166         s=abs(s)$
 167         t=scale
 168         scale=s+1
 169         v=4*a(1)
 170         scale=t
 171         return v@s
 172 }
 173 define t(x){
 174         auto s,c
 175         l=scale
 176         scale+=2
 177         s=s(x)
 178         c=c(x)
 179         scale-=2
 180         return s/c
 181 }
 182 define a2(y,x){
 183         auto a,p
 184         if(!x&&!y)y/x
 185         if(x<=0){
 186                 p=pi(scale+2)
 187                 if(y<0)p=-p
 188         }
 189         if(x==0)a=p/2
 190         else{
 191                 scale+=2
 192                 a=a(y/x)+p
 193                 scale-=2
 194         }
 195         return a@scale
 196 }
 197 define sin(x){return s(x)}
 198 define cos(x){return c(x)}
 199 define atan(x){return a(x)}
 200 define tan(x){return t(x)}
 201 define atan2(y,x){return a2(y,x)}
 202 define r2d(x){
 203         auto r,i,s
 204         s=scale
 205         scale+=5
 206         i=ibase
 207         ibase=A
 208         r=x*180/pi(scale)
 209         ibase=i
 210         scale=s
 211         return r@s
 212 }
 213 define d2r(x){
 214         auto r,i,s
 215         s=scale
 216         scale+=5
 217         i=ibase
 218         ibase=A
 219         r=x*pi(scale)/180
 220         ibase=i
 221         scale=s
 222         return r@s
 223 }
 224 define frand(p){
 225         p=abs(p)$
 226         return irand(A^p)>>p
 227 }
 228 define ifrand(i,p){return irand(abs(i)$)+frand(p)}
 229 define srand(x){
 230         if(irand(2))return -x
 231         return x
 232 }
 233 define brand(){return irand(2)}
 234 define void output(x,b){
 235         auto c
 236         c=obase
 237         obase=b
 238         x
 239         obase=c
 240 }
 241 define void hex(x){output(x,G)}
 242 define void binary(x){output(x,2)}
 243 define ubytes(x){
 244         auto p,i
 245         x=abs(x)$
 246         i=2^8
 247         for(p=1;i-1<x;p*=2){i*=i}
 248         return p
 249 }
 250 define sbytes(x){
 251         auto p,n,z
 252         z=(x<0)
 253         x=abs(x)
 254         x=x$
 255         n=ubytes(x)
 256         p=2^(n*8-1)
 257         if(x>p||(!z&&x==p))n*=2
 258         return n
 259 }
 260 define s2un(x,n){
 261         auto t,u,s
 262         x=x$
 263         if(x<0){
 264                 x=abs(x)
 265                 s=scale
 266                 scale=0
 267                 t=n*8
 268                 u=2^(t-1)
 269                 if(x==u)return x
 270                 else if(x>u)x%=u
 271                 scale=s
 272                 return 2^(t)-x
 273         }
 274         return x
 275 }
 276 define s2u(x){return s2un(x,sbytes(x))}
 277 define void plz(x){
 278         if(leading_zero())print x
 279         else{
 280                 if(x>-1&&x<1&&x!=0){
 281                         if(x<0)print"-"
 282                         print 0,abs(x)
 283                 }
 284                 else print x
 285         }
 286 }
 287 define void plznl(x){
 288         plz(x)
 289         print"\n"
 290 }
 291 define void pnlz(x){
 292         auto s,i
 293         if(leading_zero()){
 294                 if(x>-1&&x<1&&x!=0){
 295                         s=scale(x)
 296                         if(x<0)print"-"
 297                         print"."
 298                         x=abs(x)
 299                         for(i=0;i<s;++i){
 300                                 x<<=1
 301                                 print x$
 302                                 x-=x$
 303                         }
 304                         return
 305                 }
 306         }
 307         print x
 308 }
 309 define void pnlznl(x){
 310         pnlz(x)
 311         print"\n"
 312 }
 313 define void output_byte(x,i){
 314         auto j,p,y,b
 315         j=ibase
 316         ibase=A
 317         s=scale
 318         scale=0
 319         x=abs(x)$
 320         b=x/(2^(i*8))
 321         b%=256
 322         y=log(256,obase)
 323         if(b>1)p=log(b,obase)+1
 324         else p=b
 325         for(i=y-p;i>0;--i)print 0
 326         if(b)print b
 327         scale=s
 328         ibase=j
 329 }
 330 define void output_uint(x,n){
 331         auto i
 332         for(i=n-1;i>=0;--i){
 333                 output_byte(x,i)
 334                 if(i)print" "
 335                 else print"\n"
 336         }
 337 }
 338 define void hex_uint(x,n){
 339         auto o
 340         o=obase
 341         obase=G
 342         output_uint(x,n)
 343         obase=o
 344 }
 345 define void binary_uint(x,n){
 346         auto o
 347         o=obase
 348         obase=2
 349         output_uint(x,n)
 350         obase=o
 351 }
 352 define void uintn(x,n){
 353         if(scale(x)){
 354                 print"Error: ",x," is not an integer.\n"
 355                 return
 356         }
 357         if(x<0){
 358                 print"Error: ",x," is negative.\n"
 359                 return
 360         }
 361         if(x>=2^(n*8)){
 362                 print"Error: ",x," cannot fit into ",n," unsigned byte(s).\n"
 363                 return
 364         }
 365         binary_uint(x,n)
 366         hex_uint(x,n)
 367 }
 368 define void intn(x,n){
 369         auto t
 370         if(scale(x)){
 371                 print"Error: ",x," is not an integer.\n"
 372                 return
 373         }
 374         t=2^(n*8-1)
 375         if(abs(x)>=t&&(x>0||x!=-t)){
 376                 print "Error: ",x," cannot fit into ",n," signed byte(s).\n"
 377                 return
 378         }
 379         x=s2un(x,n)
 380         binary_uint(x,n)
 381         hex_uint(x,n)
 382 }
 383 define void uint8(x){uintn(x,1)}
 384 define void int8(x){intn(x,1)}
 385 define void uint16(x){uintn(x,2)}
 386 define void int16(x){intn(x,2)}
 387 define void uint32(x){uintn(x,4)}
 388 define void int32(x){intn(x,4)}
 389 define void uint64(x){uintn(x,8)}
 390 define void int64(x){intn(x,8)}
 391 define void uint(x){uintn(x,ubytes(x))}
 392 define void int(x){intn(x,sbytes(x))}
 393 define bunrev(t){
 394         auto a,s,m[]
 395         s=scale
 396         scale=0
 397         t=abs(t)$
 398         while(t!=1){
 399                 t=divmod(t,2,m[])
 400                 a*=2
 401                 a+=m[0]
 402         }
 403         scale=s
 404         return a
 405 }
 406 define band(a,b){
 407         auto s,t,m[],n[]
 408         a=abs(a)$
 409         b=abs(b)$
 410         if(b>a){
 411                 t=b
 412                 b=a
 413                 a=t
 414         }
 415         s=scale
 416         scale=0
 417         t=1
 418         while(b){
 419                 a=divmod(a,2,m[])
 420                 b=divmod(b,2,n[])
 421                 t*=2
 422                 t+=(m[0]&&n[0])
 423         }
 424         scale=s
 425         return bunrev(t)
 426 }
 427 define bor(a,b){
 428         auto s,t,m[],n[]
 429         a=abs(a)$
 430         b=abs(b)$
 431         if(b>a){
 432                 t=b
 433                 b=a
 434                 a=t
 435         }
 436         s=scale
 437         scale=0
 438         t=1
 439         while(b){
 440                 a=divmod(a,2,m[])
 441                 b=divmod(b,2,n[])
 442                 t*=2
 443                 t+=(m[0]||n[0])
 444         }
 445         while(a){
 446                 a=divmod(a,2,m[])
 447                 t*=2
 448                 t+=m[0]
 449         }
 450         scale=s
 451         return bunrev(t)
 452 }
 453 define bxor(a,b){
 454         auto s,t,m[],n[]
 455         a=abs(a)$
 456         b=abs(b)$
 457         if(b>a){
 458                 t=b
 459                 b=a
 460                 a=t
 461         }
 462         s=scale
 463         scale=0
 464         t=1
 465         while(b){
 466                 a=divmod(a,2,m[])
 467                 b=divmod(b,2,n[])
 468                 t*=2
 469                 t+=(m[0]+n[0]==1)
 470         }
 471         while(a){
 472                 a=divmod(a,2,m[])
 473                 t*=2
 474                 t+=m[0]
 475         }
 476         scale=s
 477         return bunrev(t)
 478 }
 479 define bshl(a,b){return abs(a)$*2^abs(b)$}
 480 define bshr(a,b){return (abs(a)$/2^abs(b)$)$}
 481 define bnotn(x,n){
 482         auto s,t,m[]
 483         s=scale
 484         scale=0
 485         t=2^(abs(n)$*8)
 486         x=abs(x)$%t+t
 487         t=1
 488         while(x!=1){
 489                 x=divmod(x,2,m[])
 490                 t*=2
 491                 t+=!m[0]
 492         }
 493         scale=s
 494         return bunrev(t)
 495 }
 496 define bnot8(x){return bnotn(x,1)}
 497 define bnot16(x){return bnotn(x,2)}
 498 define bnot32(x){return bnotn(x,4)}
 499 define bnot64(x){return bnotn(x,8)}
 500 define bnot(x){return bnotn(x,ubytes(x))}
 501 define brevn(x,n){
 502         auto s,t,m[]
 503         s=scale
 504         scale=0
 505         t=2^(abs(n)$*8)
 506         x=abs(x)$%t+t
 507         scale=s
 508         return bunrev(x)
 509 }
 510 define brev8(x){return brevn(x,1)}
 511 define brev16(x){return brevn(x,2)}
 512 define brev32(x){return brevn(x,4)}
 513 define brev64(x){return brevn(x,8)}
 514 define brev(x){return brevn(x,ubytes(x))}
 515 define broln(x,p,n){
 516         auto s,t,m[]
 517         s=scale
 518         scale=0
 519         n=abs(n)$*8
 520         p=abs(p)$%n
 521         t=2^n
 522         x=abs(x)$%t
 523         if(!p)return x
 524         x=divmod(x,2^(n-p),m[])
 525         x+=m[0]*2^p%t
 526         scale=s
 527         return x
 528 }
 529 define brol8(x,p){return broln(x,p,1)}
 530 define brol16(x,p){return broln(x,p,2)}
 531 define brol32(x,p){return broln(x,p,4)}
 532 define brol64(x,p){return broln(x,p,8)}
 533 define brol(x,p){return broln(x,p,ubytes(x))}
 534 define brorn(x,p,n){
 535         auto s,t,m[]
 536         s=scale
 537         scale=0
 538         n=abs(n)$*8
 539         p=abs(p)$%n
 540         t=2^n
 541         x=abs(x)$%t
 542         if(!p)return x
 543         x=divmod(x,2^p,m[])
 544         x+=m[0]*2^(n-p)%t
 545         scale=s
 546         return x
 547 }
 548 define bror8(x,p){return brorn(x,p,1)}
 549 define bror16(x,p){return brorn(x,p,2)}
 550 define bror32(x,p){return brorn(x,p,4)}
 551 define bror64(x,p){return brorn(x,p,8)}
 552 define brol(x,p){return brorn(x,p,ubytes(x))}
 553 define bmodn(x,n){
 554         auto s
 555         s=scale
 556         scale=0
 557         x=abs(x)$%2^(abs(n)$*8)
 558         scale=s
 559         return x
 560 }
 561 define bmod8(x){return bmodn(x,1)}
 562 define bmod16(x){return bmodn(x,2)}
 563 define bmod32(x){return bmodn(x,4)}
 564 define bmod64(x){return bmodn(x,8)}