1implement Pi; 2 3include "sys.m"; 4 sys: Sys; 5include "draw.m"; 6include "math.m"; 7 math: Math; 8 log: import math; 9include "daytime.m"; 10 daytime: Daytime; 11 12LBASE: con 3; # 4 13BASE: con 1000; # 10000 14 15stderr: ref Sys->FD; 16 17# spawn process for each series ? 18 19Pi: module 20{ 21 init: fn(nil: ref Draw->Context, argv: list of string); 22}; 23 24init(nil: ref Draw->Context, argv: list of string) 25{ 26 sys = load Sys Sys->PATH; 27 math = load Math Math->PATH; 28 daytime = load Daytime Daytime->PATH; 29 30 stderr = sys->fildes(2); 31 dp := 1000; 32 argv = tl argv; 33 if(argv != nil){ 34 if(tl argv != nil){ 35 picmp(hd argv, hd tl argv); 36 exit; 37 } 38 dp = int hd argv; 39 } 40 if(dp <= 0) 41 exit; 42 # t1 := daytime->now(); 43 p2 := pi2(dp+1); 44 # t2 := daytime->now(); 45 prpi(p2); 46 p1 := pi1(dp+1); 47 # t3 := daytime->now(); 48 # sys->print("%d %d\n", t2-t1, t3-t2); 49 if(p1 == nil && p2 == nil) 50 fatal("too many dp: reduce dp or source base"); 51 else if(p1 == nil) 52 p1 = p2; 53 else if(p2 == nil) 54 p2 = p1; 55 n1 := len p1; 56 n2 := len p2; 57 if(n1 != n2) 58 fatal(sys->sprint("lens differ %d %d", n1, n2)); 59 f := array[10] of { * => 0 }; 60 for(i := 0; i < n1; i++){ 61 if(p1[i] != p2[i]) 62 fatal(sys->sprint("arrays differ %d/%d: %d %d", i, n1, p1[i], p2[i])); 63 if(p1[i] < 0 || p1[i] >= BASE) 64 fatal(sys->sprint("bad array element %d: %d", i, p1[i])); 65 if(0){ 66 p := p1[i]; 67 for(j := 0; j < LBASE; j++){ 68 f[p%10]++; 69 p /= 10; 70 } 71 } 72 } 73 # prpi(p1); 74 if(0){ 75 t := 0; 76 for(i = 0; i < 10; i++){ 77 sys->print("%d %d\n", i, f[i]); 78 t += f[i]; 79 } 80 sys->print("T %d\n", t); 81 } 82} 83 84terms(dp: int, f: int, v: int): (int, int) 85{ 86 p := dp; 87 t := 0; 88 for(;;){ 89 t = 2 + int ((real p*log(real 10)+log(real v))/log(real f)); 90 if(!(t&1)) 91 t++; 92 e := int (log(real (v*(t+1)/2))/log(real 10))+1; 93 if(dp <= p-e) 94 break; 95 p += e; 96 } 97 # sys->fprint(stderr, "dp=%d p=%d f=%d v=%d terms=%d\n", dp, p, f, v, t); 98 if(t < f*f) 99 k := f*f; 100 else 101 k = t; 102 m := BASE*k; 103 if(m < 0 || m < BASE || m < k || m/BASE != k || m/k != BASE) 104 return (-1, -1); 105 return (t, p); 106} 107 108prpi(p: array of int) 109{ 110 n := len p; 111 sys->print("π ≅ "); 112 m := BASE/10; 113 sys->print("%d.%.*d", p[0]/m, LBASE-1, p[0]%m); 114 for(i := 1; i < n; i++) 115 sys->print("%.*d", LBASE, p[i]); 116 sys->print("\n"); 117} 118 119memcmp(b1: array of byte, b2: array of byte, n: int): (int, int, int) 120{ 121 for(i := 0; i < n; i++) 122 if(b1[i] != b2[i]) 123 return (i, int b1[i], int b2[i]); 124 return (-1, 0, 0); 125} 126 127picmp(f1: string, f2: string) 128{ 129 fd1 := sys->open(f1, Sys->OREAD); 130 fd2 := sys->open(f2, Sys->OREAD); 131 if(fd1 == nil || fd2 == nil) 132 fatal(sys->sprint("cannot open %s or %s", f1, f2)); 133 b1 := array[Sys->ATOMICIO] of byte; 134 b2 := array[Sys->ATOMICIO] of byte; 135 t := 0; 136 shouldexit := 0; 137 for(;;){ 138 n1 := sys->read(fd1, b1, len b1); 139 n2 := sys->read(fd2, b2, len b2); 140 if(n1 <= 0 || n2 <= 0) 141 return; 142 if(shouldexit) 143 fatal("bad picmp"); 144 if(n1 < n2) 145 (d, v1, v2) := memcmp(b1, b2, n1); 146 else 147 (d, v1, v2) = memcmp(b1, b2, n2); 148 if(d >= 0){ 149 if(v1 == '\n' || v2 == '\n') 150 shouldexit = 1; 151 else 152 fatal(sys->sprint("%s %s differ at byte %d(%c %c)", f1, f2, t+d, v1, v2)); 153 } 154 t += n1; 155 if(n1 != n2) 156 shouldexit = 1; 157 } 158} 159 160roundup(n: int, m: int): (int, int) 161{ 162 r := m*((n+m-1)/m); 163 return (r, r/m); 164} 165 166pi1(dp: int): array of int 167{ 168 fs := array[2] of { 5, 239 }; 169 vs := array[2] of { 16, 4 }; 170 ss := array[2] of { 1, -1 }; 171 # sys->fprint(stderr, "π1\n"); 172 return pi(dp, fs, vs, ss); 173} 174 175pi2(dp: int): array of int 176{ 177 fs := array[3] of { 18, 57, 239 }; 178 vs := array[3] of { 48, 32, 20 }; 179 ss := array[3] of { 1, 1, -1 }; 180 # sys->fprint(stderr, "π2\n"); 181 return pi(dp, fs, vs, ss); 182} 183 184pi3(dp: int): array of int 185{ 186 fs := array[4] of { 57, 239, 682, 12943 }; 187 vs := array[4] of { 176, 28, 48, 96 }; 188 ss := array[4] of { 1, 1, -1, 1 }; 189 # sys->fprint(stderr, "π3\n"); 190 return pi(dp, fs, vs, ss); 191} 192 193pi(dp: int, fs: array of int, vs: array of int, ss: array of int): array of int 194{ 195 k := len fs; 196 n := cn := adp := 0; 197 (dp, n) = roundup(dp, LBASE); 198 cdp := dp; 199 m := array[k] of int; 200 for(i := 0; i < k; i++){ 201 (m[i], adp) = terms(dp+1, fs[i], vs[i]); 202 if(m[i] < 0) 203 return nil; 204 if(adp > cdp) 205 cdp = adp; 206 } 207 (cdp, cn) = roundup(cdp, LBASE); 208 a := array[cn] of int; 209 p := array[cn] of int; 210 for(i = 0; i < cn; i++) 211 p[i] = 0; 212 for(i = 0; i < k; i++){ 213 series(m[i], cn, fs[i], (vs[i]*BASE)/10, ss[i], a, p); 214 # sys->fprint(stderr, "term %d done\n", i+1); 215 } 216 return p[0: n]; 217} 218 219series(m: int, n: int, f: int, v: int, s: int, a: array of int, p: array of int) 220{ 221 i, j, k, q, r, r1, r2, n0: int; 222 223 v *= f; 224 f *= f; 225 for(j = 0; j < n; j++) 226 a[j] = 0; 227 a[0] = v; 228 229 if(s == 1) 230 series1(m, n, f, v, a, p); 231 else 232 series2(m, n, f, v, a, p); 233 return; 234 235 # following code now split 236 n0 = 0; # reaches n when very close to m so no check needed 237 for(i = 1; i <= m; i += 2){ 238 r1 = r2 = 0; 239 for(j = n0; j < n; j++){ 240 v = a[j]+r1; 241 q = v/f; 242 r1 = (v-q*f)*BASE; 243 a[j] = q; 244 v = q+r2; 245 q = v/i; 246 r2 = (v-q*i)*BASE; 247 for(k = j; q > 0; k--){ 248 r = p[k]+s*q; 249 if(r >= BASE){ 250 p[k] = r-BASE; 251 q = 1; 252 } 253 else if(r < 0){ 254 p[k] = r+BASE; 255 q = 1; 256 } 257 else{ 258 p[k] = r; 259 q = 0; 260 } 261 } 262 } 263 for(j = n0; j < n; j++){ 264 if(a[j] == 0) 265 n0++; 266 else 267 break; 268 } 269 s = -s; 270 } 271} 272 273series1(m: int, n: int, f: int, v: int, a: array of int, p: array of int) 274{ 275 i, j, k, q, r, r1, r2, n0: int; 276 277 n0 = 0; 278 for(i = 1; i <= m; i += 2){ 279 r1 = r2 = 0; 280 for(j = n0; j < n; j++){ 281 v = a[j]+r1; 282 q = v/f; 283 r1 = (v-q*f)*BASE; 284 a[j] = q; 285 v = q+r2; 286 q = v/i; 287 r2 = (v-q*i)*BASE; 288 for(k = j; q > 0; k--){ 289 r = p[k]+q; 290 if(r >= BASE){ 291 p[k] = r-BASE; 292 q = 1; 293 } 294 else{ 295 p[k] = r; 296 q = 0; 297 } 298 } 299 } 300 for(j = n0; j < n; j++){ 301 if(a[j] == 0) 302 n0++; 303 else 304 break; 305 } 306 i += 2; 307 r1 = r2 = 0; 308 for(j = n0; j < n; j++){ 309 v = a[j]+r1; 310 q = v/f; 311 r1 = (v-q*f)*BASE; 312 a[j] = q; 313 v = q+r2; 314 q = v/i; 315 r2 = (v-q*i)*BASE; 316 for(k = j; q > 0; k--){ 317 r = p[k]-q; 318 if(r < 0){ 319 p[k] = r+BASE; 320 q = 1; 321 } 322 else{ 323 p[k] = r; 324 q = 0; 325 } 326 } 327 } 328 for(j = n0; j < n; j++){ 329 if(a[j] == 0) 330 n0++; 331 else 332 break; 333 } 334 } 335} 336 337series2(m: int, n: int, f: int, v: int, a: array of int, p: array of int) 338{ 339 i, j, k, q, r, r1, r2, n0: int; 340 341 n0 = 0; 342 for(i = 1; i <= m; i += 2){ 343 r1 = r2 = 0; 344 for(j = n0; j < n; j++){ 345 v = a[j]+r1; 346 q = v/f; 347 r1 = (v-q*f)*BASE; 348 a[j] = q; 349 v = q+r2; 350 q = v/i; 351 r2 = (v-q*i)*BASE; 352 for(k = j; q > 0; k--){ 353 r = p[k]-q; 354 if(r < 0){ 355 p[k] = r+BASE; 356 q = 1; 357 } 358 else{ 359 p[k] = r; 360 q = 0; 361 } 362 } 363 } 364 for(j = n0; j < n; j++){ 365 if(a[j] == 0) 366 n0++; 367 else 368 break; 369 } 370 i += 2; 371 r1 = r2 = 0; 372 for(j = n0; j < n; j++){ 373 v = a[j]+r1; 374 q = v/f; 375 r1 = (v-q*f)*BASE; 376 a[j] = q; 377 v = q+r2; 378 q = v/i; 379 r2 = (v-q*i)*BASE; 380 for(k = j; q > 0; k--){ 381 r = p[k]+q; 382 if(r >= BASE){ 383 p[k] = r-BASE; 384 q = 1; 385 } 386 else{ 387 p[k] = r; 388 q = 0; 389 } 390 } 391 } 392 for(j = n0; j < n; j++){ 393 if(a[j] == 0) 394 n0++; 395 else 396 break; 397 } 398 } 399} 400 401fatal(e: string) 402{ 403 sys->print("%s\n", e); 404 exit; 405} 406