contrib/bc/gen/lib2.bc

   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 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         n=ubytes(x)
 255         p=2^(n*8-1)
 256         if(x>p||(!z&&x==p))n*=2
 257         return n
 258 }
 259 define s2un(x,n){
 260         auto t,u,s
 261         x=x$
 262         if(x<0){
 263                 x=abs(x)
 264                 s=scale
 265                 scale=0
 266                 t=n*8
 267                 u=2^(t-1)
 268                 if(x==u)return x
 269                 else if(x>u)x%=u
 270                 scale=s
 271                 return 2^(t)-x
 272         }
 273         return x
 274 }
 275 define s2u(x){return s2un(x,sbytes(x))}
 276 define void plz(x){
 277         if(leading_zero())print x
 278         else{
 279                 if(x>-1&&x<1&&x!=0){
 280                         if(x<0)print"-"
 281                         print 0,abs(x)
 282                 }
 283                 else print x
 284         }
 285 }
 286 define void plznl(x){
 287         plz(x)
 288         print"\n"
 289 }
 290 define void pnlz(x){
 291         auto s,i
 292         if(leading_zero()){
 293                 if(x>-1&&x<1&&x!=0){
 294                         s=scale(x)
 295                         if(x<0)print"-"
 296                         print"."
 297                         x=abs(x)
 298                         for(i=0;i<s;++i){
 299                                 x<<=1
 300                                 print x$
 301                                 x-=x$
 302                         }
 303                         return
 304                 }
 305         }
 306         print x
 307 }
 308 define void pnlznl(x){
 309         pnlz(x)
 310         print"\n"
 311 }
 312 define void output_byte(x,i){
 313         auto j,p,y,b,s
 314         s=scale
 315         scale=0
 316         x=abs(x)$
 317         b=x/(2^(i*8))
 318         j=2^8
 319         b%=j
 320         y=log(j,obase)
 321         if(b>1)p=log(b,obase)+1
 322         else p=b
 323         for(i=y-p;i>0;--i)print 0
 324         if(b)print b
 325         scale=s
 326 }
 327 define void output_uint(x,n){
 328         auto i
 329         for(i=n-1;i>=0;--i){
 330                 output_byte(x,i)
 331                 if(i)print" "
 332                 else print"\n"
 333         }
 334 }
 335 define void hex_uint(x,n){
 336         auto o
 337         o=obase
 338         obase=G
 339         output_uint(x,n)
 340         obase=o
 341 }
 342 define void binary_uint(x,n){
 343         auto o
 344         o=obase
 345         obase=2
 346         output_uint(x,n)
 347         obase=o
 348 }
 349 define void uintn(x,n){
 350         if(scale(x)){
 351                 print"Error: ",x," is not an integer.\n"
 352                 return
 353         }
 354         if(x<0){
 355                 print"Error: ",x," is negative.\n"
 356                 return
 357         }
 358         if(x>=2^(n*8)){
 359                 print"Error: ",x," cannot fit into ",n," unsigned byte(s).\n"
 360                 return
 361         }
 362         binary_uint(x,n)
 363         hex_uint(x,n)
 364 }
 365 define void intn(x,n){
 366         auto t
 367         if(scale(x)){
 368                 print"Error: ",x," is not an integer.\n"
 369                 return
 370         }
 371         t=2^(n*8-1)
 372         if(abs(x)>=t&&(x>0||x!=-t)){
 373                 print "Error: ",x," cannot fit into ",n," signed byte(s).\n"
 374                 return
 375         }
 376         x=s2un(x,n)
 377         binary_uint(x,n)
 378         hex_uint(x,n)
 379 }
 380 define void uint8(x){uintn(x,1)}
 381 define void int8(x){intn(x,1)}
 382 define void uint16(x){uintn(x,2)}
 383 define void int16(x){intn(x,2)}
 384 define void uint32(x){uintn(x,4)}
 385 define void int32(x){intn(x,4)}
 386 define void uint64(x){uintn(x,8)}
 387 define void int64(x){intn(x,8)}
 388 define void uint(x){uintn(x,ubytes(x))}
 389 define void int(x){intn(x,sbytes(x))}
 390 define bunrev(t){
 391         auto a,s,m[]
 392         s=scale
 393         scale=0
 394         t=abs(t)$
 395         while(t!=1){
 396                 t=divmod(t,2,m[])
 397                 a*=2
 398                 a+=m[0]
 399         }
 400         scale=s
 401         return a
 402 }
 403 define band(a,b){
 404         auto s,t,m[],n[]
 405         a=abs(a)$
 406         b=abs(b)$
 407         if(b>a){
 408                 t=b
 409                 b=a
 410                 a=t
 411         }
 412         s=scale
 413         scale=0
 414         t=1
 415         while(b){
 416                 a=divmod(a,2,m[])
 417                 b=divmod(b,2,n[])
 418                 t*=2
 419                 t+=(m[0]&&n[0])
 420         }
 421         scale=s
 422         return bunrev(t)
 423 }
 424 define bor(a,b){
 425         auto s,t,m[],n[]
 426         a=abs(a)$
 427         b=abs(b)$
 428         if(b>a){
 429                 t=b
 430                 b=a
 431                 a=t
 432         }
 433         s=scale
 434         scale=0
 435         t=1
 436         while(b){
 437                 a=divmod(a,2,m[])
 438                 b=divmod(b,2,n[])
 439                 t*=2
 440                 t+=(m[0]||n[0])
 441         }
 442         while(a){
 443                 a=divmod(a,2,m[])
 444                 t*=2
 445                 t+=m[0]
 446         }
 447         scale=s
 448         return bunrev(t)
 449 }
 450 define bxor(a,b){
 451         auto s,t,m[],n[]
 452         a=abs(a)$
 453         b=abs(b)$
 454         if(b>a){
 455                 t=b
 456                 b=a
 457                 a=t
 458         }
 459         s=scale
 460         scale=0
 461         t=1
 462         while(b){
 463                 a=divmod(a,2,m[])
 464                 b=divmod(b,2,n[])
 465                 t*=2
 466                 t+=(m[0]+n[0]==1)
 467         }
 468         while(a){
 469                 a=divmod(a,2,m[])
 470                 t*=2
 471                 t+=m[0]
 472         }
 473         scale=s
 474         return bunrev(t)
 475 }
 476 define bshl(a,b){return abs(a)$*2^abs(b)$}
 477 define bshr(a,b){return (abs(a)$/2^abs(b)$)$}
 478 define bnotn(x,n){
 479         auto s,t,m[]
 480         s=scale
 481         scale=0
 482         t=2^(abs(n)$*8)
 483         x=abs(x)$%t+t
 484         t=1
 485         while(x!=1){
 486                 x=divmod(x,2,m[])
 487                 t*=2
 488                 t+=!m[0]
 489         }
 490         scale=s
 491         return bunrev(t)
 492 }
 493 define bnot8(x){return bnotn(x,1)}
 494 define bnot16(x){return bnotn(x,2)}
 495 define bnot32(x){return bnotn(x,4)}
 496 define bnot64(x){return bnotn(x,8)}
 497 define bnot(x){return bnotn(x,ubytes(x))}
 498 define brevn(x,n){
 499         auto s,t,m[]
 500         s=scale
 501         scale=0
 502         t=2^(abs(n)$*8)
 503         x=abs(x)$%t+t
 504         scale=s
 505         return bunrev(x)
 506 }
 507 define brev8(x){return brevn(x,1)}
 508 define brev16(x){return brevn(x,2)}
 509 define brev32(x){return brevn(x,4)}
 510 define brev64(x){return brevn(x,8)}
 511 define brev(x){return brevn(x,ubytes(x))}
 512 define broln(x,p,n){
 513         auto s,t,m[]
 514         s=scale
 515         scale=0
 516         n=abs(n)$*8
 517         p=abs(p)$%n
 518         t=2^n
 519         x=abs(x)$%t
 520         if(!p)return x
 521         x=divmod(x,2^(n-p),m[])
 522         x+=m[0]*2^p%t
 523         scale=s
 524         return x
 525 }
 526 define brol8(x,p){return broln(x,p,1)}
 527 define brol16(x,p){return broln(x,p,2)}
 528 define brol32(x,p){return broln(x,p,4)}
 529 define brol64(x,p){return broln(x,p,8)}
 530 define brol(x,p){return broln(x,p,ubytes(x))}
 531 define brorn(x,p,n){
 532         auto s,t,m[]
 533         s=scale
 534         scale=0
 535         n=abs(n)$*8
 536         p=abs(p)$%n
 537         t=2^n
 538         x=abs(x)$%t
 539         if(!p)return x
 540         x=divmod(x,2^p,m[])
 541         x+=m[0]*2^(n-p)%t
 542         scale=s
 543         return x
 544 }
 545 define bror8(x,p){return brorn(x,p,1)}
 546 define bror16(x,p){return brorn(x,p,2)}
 547 define bror32(x,p){return brorn(x,p,4)}
 548 define bror64(x,p){return brorn(x,p,8)}
 549 define brol(x,p){return brorn(x,p,ubytes(x))}
 550 define bmodn(x,n){
 551         auto s
 552         s=scale
 553         scale=0
 554         x=abs(x)$%2^(abs(n)$*8)
 555         scale=s
 556         return x
 557 }
 558 define bmod8(x){return bmodn(x,1)}
 559 define bmod16(x){return bmodn(x,2)}
 560 define bmod32(x){return bmodn(x,4)}
 561 define bmod64(x){return bmodn(x,8)}