File: common\fractalb.c
1 /* -----------------------------------------------------------------
2
3 This file contains the "big number" high precision versions of the
4 fractal routines.
5
6 -------------------------------------------------------------------- */
7
8
9 #include <limits.h>
10 #include <string.h>
11 #ifdef __TURBOC__
13 #elif !defined(__386BSD__)
14 #include <malloc.h>
15 #endif
16 /* see Fractint.c for a description of the "include" hierarchy */
17 #include "port.h"
18 #include "prototyp.h"
19 #include "helpdefs.h"
20 #include "fractype.h"
21
22
23 int bf_math = 0;
24
25 #if (_MSC_VER >= 700)
26 #pragma code_seg ("bigsetup_text") /* place following in an overlay */
27 #endif
28
29 #ifdef DEBUG
30
31 /**********************************************************************/
32 void show_var_bn(char *s, bn_t n)
33 {
34 char msg[200];
35 strcpy(msg,s);
36 strcat(msg," ");
37 bntostr(msg+strlen(s),40,n);
38 msg[79] = 0;
39 stopmsg(0,(char far *)msg);
40 }
41
42 void showcornersdbl(char *s)
43 {
44 char msg[400];
45 sprintf(msg,"%s\n\
46 xxmin= %.20f xxmax= %.20f\n\
47 yymin= %.20f yymax= %.20f\n\
48 xx3rd= %.20f yy3rd= %.20f\n\
49 delxx= %.20Lf delyy= %.20Lf\n\
50 delx2= %.20Lf dely2= %.20Lf",s,xxmin,xxmax,yymin,yymax,xx3rd,yy3rd,
51 delxx, delyy,delxx2, delyy2);
52 if(stopmsg(0,msg)==-1)
53 goodbye();
54 }
55
56 /* show floating point and bignumber corners */
57 void showcorners(char *s)
58 {
59 int dec=20;
60 char msg[100],msg1[100],msg3[100];
61 bntostr(msg,dec,bnxmin);
62 sprintf(msg1,"bnxmin=%s\nxxmin= %.20f\n\n",msg,xxmin);
63 strcpy(msg3,s);
64 strcat(msg3,"\n");
65 strcat(msg3,msg1);
66 bntostr(msg,dec,bnxmax);
67 sprintf(msg1,"bnxmax=%s\nxxmax= %.20f\n\n",msg,xxmax);
68 strcat(msg3,msg1);
69 bntostr(msg,dec,bnymin);
70 sprintf(msg1,"bnymin=%s\nyymin= %.20f\n\n",msg,yymin);
71 strcat(msg3,msg1);
72 bntostr(msg,dec,bnymax);
73 sprintf(msg1,"bnymax=%s\nyymax= %.20f\n\n",msg,yymax);
74 strcat(msg3,msg1);
75 bntostr(msg,dec,bnx3rd);
76 sprintf(msg1,"bnx3rd=%s\nxx3rd= %.20f\n\n",msg,xx3rd);
77 strcat(msg3,msg1);
78 bntostr(msg,dec,bny3rd);
79 sprintf(msg1,"bny3rd=%s\nyy3rd= %.20f\n\n",msg,yy3rd);
80 strcat(msg3,msg1);
81 if(stopmsg(0,msg3)==-1)
82 goodbye();
83 }
84
85 /* show globals */
86 void showbfglobals(char *s)
87 {
88 char msg[300];
89 sprintf(msg, "%s\n\
90 bnstep=%d bnlength=%d intlength=%d rlength=%d padding=%d\n\
91 shiftfactor=%d decimals=%d bflength=%d rbflength=%d \n\
92 bfdecimals=%d ",
93 s, bnstep, bnlength, intlength, rlength, padding,
94 shiftfactor, decimals, bflength, rbflength,
95 bfdecimals);
96 if(stopmsg(0,msg)==-1)
97 goodbye();
98 }
99
100 void showcornersbf(char *s)
101 {
102 int dec=decimals;
103 char msg[100],msg1[100],msg3[600];
104 if(dec > 20) dec = 20;
105 bftostr(msg,dec,bfxmin);
106 sprintf(msg1,"bfxmin=%s\nxxmin= %.20f decimals %d bflength %d\n\n",
107 msg,xxmin,decimals,bflength);
108 strcpy(msg3,s);
109 strcat(msg3,"\n");
110 strcat(msg3,msg1);
111 bftostr(msg,dec,bfxmax);
112 sprintf(msg1,"bfxmax=%s\nxxmax= %.20f\n\n",msg,xxmax);
113 strcat(msg3,msg1);
114 bftostr(msg,dec,bfymin);
115 sprintf(msg1,"bfymin=%s\nyymin= %.20f\n\n",msg,yymin);
116 strcat(msg3,msg1);
117 bftostr(msg,dec,bfymax);
118 sprintf(msg1,"bfymax=%s\nyymax= %.20f\n\n",msg,yymax);
119 strcat(msg3,msg1);
120 bftostr(msg,dec,bfx3rd);
121 sprintf(msg1,"bfx3rd=%s\nxx3rd= %.20f\n\n",msg,xx3rd);
122 strcat(msg3,msg1);
123 bftostr(msg,dec,bfy3rd);
124 sprintf(msg1,"bfy3rd=%s\nyy3rd= %.20f\n\n",msg,yy3rd);
125 strcat(msg3,msg1);
126 if(stopmsg(0,msg3)==-1)
127 goodbye();
128 }
129
130 void showcornersbfs(char *s)
131 {
132 int dec=20;
133 char msg[100],msg1[100],msg3[500];
134 bftostr(msg,dec,bfsxmin);
135 sprintf(msg1,"bfsxmin=%s\nxxmin= %.20f\n\n",msg,xxmin);
136 strcpy(msg3,s);
137 strcat(msg3,"\n");
138 strcat(msg3,msg1);
139 bftostr(msg,dec,bfsxmax);
140 sprintf(msg1,"bfsxmax=%s\nxxmax= %.20f\n\n",msg,xxmax);
141 strcat(msg3,msg1);
142 bftostr(msg,dec,bfsymin);
143 sprintf(msg1,"bfsymin=%s\nyymin= %.20f\n\n",msg,yymin);
144 strcat(msg3,msg1);
145 bftostr(msg,dec,bfsymax);
146 sprintf(msg1,"bfsymax=%s\nyymax= %.20f\n\n",msg,yymax);
147 strcat(msg3,msg1);
148 bftostr(msg,dec,bfsx3rd);
149 sprintf(msg1,"bfsx3rd=%s\nxx3rd= %.20f\n\n",msg,xx3rd);
150 strcat(msg3,msg1);
151 bftostr(msg,dec,bfsy3rd);
152 sprintf(msg1,"bfsy3rd=%s\nyy3rd= %.20f\n\n",msg,yy3rd);
153 strcat(msg3,msg1);
154 if(stopmsg(0,msg3)==-1)
155 goodbye();
156 }
157
158 void show_two_bf(char *s1,bf_t t1,char *s2, bf_t t2, int digits)
159 {
160 char msg1[200],msg2[200], msg3[400];
161 bftostr_e(msg1,digits,t1);
162 bftostr_e(msg2,digits,t2);
163 sprintf(msg3,"\n%s->%s\n%s->%s",s1,msg1,s2,msg2);
164 if(stopmsg(0,msg3)==-1)
165 goodbye();
166 }
167
168 void show_three_bf(char *s1,bf_t t1,char *s2, bf_t t2, char *s3, bf_t t3, int digits)
169 {
170 char msg1[200],msg2[200], msg3[200], msg4[600];
171 bftostr_e(msg1,digits,t1);
172 bftostr_e(msg2,digits,t2);
173 bftostr_e(msg3,digits,t3);
174 sprintf(msg4,"\n%s->%s\n%s->%s\n%s->%s",s1,msg1,s2,msg2,s3,msg3);
175 if(stopmsg(0,msg4)==-1)
176 goodbye();
177 }
178
179 /* for aspect ratio debugging */
180 void showaspect(char *s)
181 {
182 bf_t bt1,bt2,aspect;
183 char msg[100],str[100];
184 int saved; saved = save_stack();
185 bt1 = alloc_stack(rbflength+2);
186 bt2 = alloc_stack(rbflength+2);
187 aspect = alloc_stack(rbflength+2);
188 sub_bf(bt1,bfxmax,bfxmin);
189 sub_bf(bt2,bfymax,bfymin);
190 div_bf(aspect,bt2,bt1);
191 bftostr(str,10,aspect);
192 sprintf(msg,"aspect %s\nfloat %13.10f\nbf %s\n\n",
193 s,
194 (yymax-yymin)/(xxmax-xxmin),
195 str);
196 if(stopmsg(0,msg)==-1)
197 goodbye();
198 restore_stack(saved);
199 }
200
201 /* compare a double and bignumber */
202 void comparevalues(char *s, LDBL x, bn_t bnx)
203 {
204 int dec=40;
205 char msg[100],msg1[100];
206 bntostr(msg,dec,bnx);
207 sprintf(msg1,"%s\nbignum=%s\ndouble=%.20Lf\n\n",s,msg,x);
208 if(stopmsg(0,msg1)==-1)
209 goodbye();
210 }
211 /* compare a double and bignumber */
212 void comparevaluesbf(char *s, LDBL x, bf_t bfx)
213 {
214 int dec=40;
215 char msg[300],msg1[300];
216 bftostr_e(msg,dec,bfx);
217 sprintf(msg1,"%s\nbignum=%s\ndouble=%.20Lf\n\n",s,msg,x);
218 if(stopmsg(0,msg1)==-1)
219 goodbye();
220 }
221
222 /**********************************************************************/
223 void show_var_bf(char *s, bf_t n)
224 {
225 char msg[200];
226 strcpy(msg,s);
227 strcat(msg," ");
228 bftostr_e(msg+strlen(s),40,n);
229 msg[79] = 0;
230 if(stopmsg(0,msg)==-1)
231 goodbye();
232 }
233 234 #endif
235
236 void bfcornerstofloat(void)
237 {
238 int i;
239 if(bf_math)
240 {
241 xxmin = (double)bftofloat(bfxmin);
242 yymin = (double)bftofloat(bfymin);
243 xxmax = (double)bftofloat(bfxmax);
244 yymax = (double)bftofloat(bfymax);
245 xx3rd = (double)bftofloat(bfx3rd);
246 yy3rd = (double)bftofloat(bfy3rd);
247 }
248 for(i=0;i<MAXPARAMS;i++)
249 if(typehasparm(fractype,i,NULL))
250 param[i] = (double)bftofloat(bfparms[i]);
251 }
252
253 #if (_MSC_VER >= 700)
254 #pragma code_seg ( ) /* back to normal segment */
255 #endif
256
257 /* -------------------------------------------------------------------- */
258 /* Bignumber Bailout Routines */
259 /* -------------------------------------------------------------------- */
260
261 /* mandel_bntoint() can only be used for intlength of 1 */
262 #define mandel_bntoint(n) (*(n + bnlength - 1)) /* assumes intlength of 1 */
263
264 /* Note: */
265 /* No need to set magnitude */
266 /* as color schemes that need it calculate it later. */
267
268 int near bnMODbailout()
269 {
270 long longmagnitude;
271
272 square_bn(bntmpsqrx, bnnew.x);
273 square_bn(bntmpsqry, bnnew.y);
274 add_bn(bntmp, bntmpsqrx+shiftfactor, bntmpsqry+shiftfactor);
275
276 longmagnitude = bntoint(bntmp); /* works with any fractal type */
277 if (longmagnitude >= (long)rqlim)
278 return 1;
279 copy_bn(bnold.x, bnnew.x);
280 copy_bn(bnold.y, bnnew.y);
281 return 0;
282 }
283
284 int near bnREALbailout()
285 {
286 long longtempsqrx;
287
288 square_bn(bntmpsqrx, bnnew.x);
289 square_bn(bntmpsqry, bnnew.y);
290 longtempsqrx = bntoint(bntmpsqrx+shiftfactor);
291 if (longtempsqrx >= (long)rqlim)
292 return 1;
293 copy_bn(bnold.x, bnnew.x);
294 copy_bn(bnold.y, bnnew.y);
295 return 0;
296 }
297
298
299 int near bnIMAGbailout()
300 {
301 long longtempsqry;
302
303 square_bn(bntmpsqrx, bnnew.x);
304 square_bn(bntmpsqry, bnnew.y);
305 longtempsqry = bntoint(bntmpsqry+shiftfactor);
306 if (longtempsqry >= (long)rqlim)
307 return 1;
308 copy_bn(bnold.x, bnnew.x);
309 copy_bn(bnold.y, bnnew.y);
310 return(0);
311 }
312
313 int near bnORbailout()
314 {
315 long longtempsqrx, longtempsqry;
316
317 square_bn(bntmpsqrx, bnnew.x);
318 square_bn(bntmpsqry, bnnew.y);
319 longtempsqrx = bntoint(bntmpsqrx+shiftfactor);
320 longtempsqry = bntoint(bntmpsqry+shiftfactor);
321 if(longtempsqrx >= (long)rqlim || longtempsqry >= (long)rqlim)
322 return 1;
323 copy_bn(bnold.x, bnnew.x);
324 copy_bn(bnold.y, bnnew.y);
325 return(0);
326 }
327
328 int near bnANDbailout()
329 {
330 long longtempsqrx, longtempsqry;
331
332 square_bn(bntmpsqrx, bnnew.x);
333 square_bn(bntmpsqry, bnnew.y);
334 longtempsqrx = bntoint(bntmpsqrx+shiftfactor);
335 longtempsqry = bntoint(bntmpsqry+shiftfactor);
336 if(longtempsqrx >= (long)rqlim && longtempsqry >= (long)rqlim)
337 return 1;
338 copy_bn(bnold.x, bnnew.x);
339 copy_bn(bnold.y, bnnew.y);
340 return(0);
341 }
342
343 int near bnMANHbailout()
344 {
345 long longtempmag;
346
347 square_bn(bntmpsqrx, bnnew.x);
348 square_bn(bntmpsqry, bnnew.y);
349 /* note: in next five lines, bnold is just used as a temporary variable */
350 abs_bn(bnold.x,bnnew.x);
351 abs_bn(bnold.y,bnnew.y);
352 add_bn(bntmp, bnold.x, bnold.y);
353 square_bn(bnold.x, bntmp);
354 longtempmag = bntoint(bnold.x+shiftfactor);
355 if(longtempmag >= (long)rqlim)
356 return 1;
357 copy_bn(bnold.x, bnnew.x);
358 copy_bn(bnold.y, bnnew.y);
359 return(0);
360 }
361
362 int near bnMANRbailout()
363 {
364 long longtempmag;
365
366 square_bn(bntmpsqrx, bnnew.x);
367 square_bn(bntmpsqry, bnnew.y);
368 add_bn(bntmp, bnnew.x, bnnew.y); /* don't need abs since we square it next */
369 /* note: in next two lines, bnold is just used as a temporary variable */
370 square_bn(bnold.x, bntmp);
371 longtempmag = bntoint(bnold.x+shiftfactor);
372 if(longtempmag >= (long)rqlim)
373 return 1;
374 copy_bn(bnold.x, bnnew.x);
375 copy_bn(bnold.y, bnnew.y);
376 return(0);
377 }
378
379 int near bfMODbailout()
380 {
381 long longmagnitude;
382
383 square_bf(bftmpsqrx, bfnew.x);
384 square_bf(bftmpsqry, bfnew.y);
385 add_bf(bftmp, bftmpsqrx, bftmpsqry);
386
387 longmagnitude = bftoint(bftmp);
388 if (longmagnitude >= (long)rqlim)
389 return 1;
390 copy_bf(bfold.x, bfnew.x);
391 copy_bf(bfold.y, bfnew.y);
392 return 0;
393 }
394
395 int near bfREALbailout()
396 {
397 long longtempsqrx;
398
399 square_bf(bftmpsqrx, bfnew.x);
400 square_bf(bftmpsqry, bfnew.y);
401 longtempsqrx = bftoint(bftmpsqrx);
402 if (longtempsqrx >= (long)rqlim)
403 return 1;
404 copy_bf(bfold.x, bfnew.x);
405 copy_bf(bfold.y, bfnew.y);
406 return 0;
407 }
408
409
410 int near bfIMAGbailout()
411 {
412 long longtempsqry;
413
414 square_bf(bftmpsqrx, bfnew.x);
415 square_bf(bftmpsqry, bfnew.y);
416 longtempsqry = bftoint(bftmpsqry);
417 if (longtempsqry >= (long)rqlim)
418 return 1;
419 copy_bf(bfold.x, bfnew.x);
420 copy_bf(bfold.y, bfnew.y);
421 return(0);
422 }
423
424 int near bfORbailout()
425 {
426 long longtempsqrx, longtempsqry;
427
428 square_bf(bftmpsqrx, bfnew.x);
429 square_bf(bftmpsqry, bfnew.y);
430 longtempsqrx = bftoint(bftmpsqrx);
431 longtempsqry = bftoint(bftmpsqry);
432 if(longtempsqrx >= (long)rqlim || longtempsqry >= (long)rqlim)
433 return 1;
434 copy_bf(bfold.x, bfnew.x);
435 copy_bf(bfold.y, bfnew.y);
436 return(0);
437 }
438
439 int near bfANDbailout()
440 {
441 long longtempsqrx, longtempsqry;
442
443 square_bf(bftmpsqrx, bfnew.x);
444 square_bf(bftmpsqry, bfnew.y);
445 longtempsqrx = bftoint(bftmpsqrx);
446 longtempsqry = bftoint(bftmpsqry);
447 if(longtempsqrx >= (long)rqlim && longtempsqry >= (long)rqlim)
448 return 1;
449 copy_bf(bfold.x, bfnew.x);
450 copy_bf(bfold.y, bfnew.y);
451 return(0);
452 }
453
454 int near bfMANHbailout()
455 {
456 long longtempmag;
457
458 square_bf(bftmpsqrx, bfnew.x);
459 square_bf(bftmpsqry, bfnew.y);
460 /* note: in next five lines, bfold is just used as a temporary variable */
461 abs_bf(bfold.x,bfnew.x);
462 abs_bf(bfold.y,bfnew.y);
463 add_bf(bftmp, bfold.x, bfold.y);
464 square_bf(bfold.x, bftmp);
465 longtempmag = bftoint(bfold.x);
466 if(longtempmag >= (long)rqlim)
467 return 1;
468 copy_bf(bfold.x, bfnew.x);
469 copy_bf(bfold.y, bfnew.y);
470 return(0);
471 }
472
473 int near bfMANRbailout()
474 {
475 long longtempmag;
476
477 square_bf(bftmpsqrx, bfnew.x);
478 square_bf(bftmpsqry, bfnew.y);
479 add_bf(bftmp, bfnew.x, bfnew.y); /* don't need abs since we square it next */
480 /* note: in next two lines, bfold is just used as a temporary variable */
481 square_bf(bfold.x, bftmp);
482 longtempmag = bftoint(bfold.x);
483 if(longtempmag >= (long)rqlim)
484 return 1;
485 copy_bf(bfold.x, bfnew.x);
486 copy_bf(bfold.y, bfnew.y);
487 return(0);
488 }
489
490 #if (_MSC_VER >= 700)
491 #pragma code_seg ("bigsetup_text") /* place following in an overlay */
492 #endif
493
494 int MandelbnSetup()
495 {
496 /* this should be set up dynamically based on corners */
497 bn_t bntemp1, bntemp2;
498 int saved; saved = save_stack();
499 bntemp1 = alloc_stack(bnlength);
500 bntemp2 = alloc_stack(bnlength);
501
502 bftobn(bnxmin,bfxmin);
503 bftobn(bnxmax,bfxmax);
504 bftobn(bnymin,bfymin);
505 bftobn(bnymax,bfymax);
506 bftobn(bnx3rd,bfx3rd);
507 bftobn(bny3rd,bfy3rd);
508
509 bf_math = BIGNUM;
510
511 /* bnxdel = (bnxmax - bnx3rd)/(xdots-1) */
512 sub_bn(bnxdel, bnxmax, bnx3rd);
513 div_a_bn_int(bnxdel, (U16)(xdots - 1));
514
515 /* bnydel = (bnymax - bny3rd)/(ydots-1) */
516 sub_bn(bnydel, bnymax, bny3rd);
517 div_a_bn_int(bnydel, (U16)(ydots - 1));
518
519 /* bnxdel2 = (bnx3rd - bnxmin)/(ydots-1) */
520 sub_bn(bnxdel2, bnx3rd, bnxmin);
521 div_a_bn_int(bnxdel2, (U16)(ydots - 1));
522
523 /* bnydel2 = (bny3rd - bnymin)/(xdots-1) */
524 sub_bn(bnydel2, bny3rd, bnymin);
525 div_a_bn_int(bnydel2, (U16)(xdots - 1));
526
527 abs_bn(bnclosenuff,bnxdel);
528 if(cmp_bn(abs_bn(bntemp1,bnxdel2),bnclosenuff) > 0)
529 copy_bn(bnclosenuff,bntemp1);
530 if(cmp_bn(abs_bn(bntemp1,bnydel),abs_bn(bntemp2,bnydel2)) > 0)
531 {
532 if(cmp_bn(bntemp1,bnclosenuff) > 0)
533 copy_bn(bnclosenuff,bntemp1);
534 }
535 else if(cmp_bn(bntemp2,bnclosenuff) > 0)
536 copy_bn(bnclosenuff,bntemp2);
537 {
538 int t;
539 t = abs(periodicitycheck);
540 while(t--)
541 half_a_bn(bnclosenuff);
542 }
543
544 c_exp = (int)param[2];
545 switch (fractype)
546 {
547 case JULIAFP:
548 bftobn(bnparm.x, bfparms[0]);
549 bftobn(bnparm.y, bfparms[1]);
550 break;
551 case FPMANDELZPOWER:
552 init_big_pi();
553 if((double)c_exp == param[2] && (c_exp & 1)) /* odd exponents */
554 symmetry = XYAXIS_NOPARM;
555 if(param[3] != 0)
556 symmetry = NOSYM;
557 break;
558 case FPJULIAZPOWER:
559 init_big_pi();
560 bftobn(bnparm.x, bfparms[0]);
561 bftobn(bnparm.y, bfparms[1]);
562 if((c_exp & 1) || param[3] != 0.0 || (double)c_exp != param[2] )
563 symmetry = NOSYM;
564 break;
565 }
566
567 /* at the present time, parameters are kept in float, but want to keep
568 the arbitrary precision logic intact. The next two lines, if used,
569 would disguise and breaking of the arbitrary precision logic */
570 /*
571 floattobn(bnparm.x,param[0]);
572 floattobn(bnparm.y,param[1]);
573 */
574 restore_stack(saved);
575 return (1);
576 }
577
578 int MandelbfSetup()
579 {
580 /* this should be set up dynamically based on corners */
581 bf_t bftemp1, bftemp2;
582 int saved; saved = save_stack();
583 bftemp1 = alloc_stack(bflength+2);
584 bftemp2 = alloc_stack(bflength+2);
585
586 bf_math = BIGFLT;
587
588 /* bfxdel = (bfxmax - bfx3rd)/(xdots-1) */
589 sub_bf(bfxdel, bfxmax, bfx3rd);
590 div_a_bf_int(bfxdel, (U16)(xdots - 1));
591
592 /* bfydel = (bfymax - bfy3rd)/(ydots-1) */
593 sub_bf(bfydel, bfymax, bfy3rd);
594 div_a_bf_int(bfydel, (U16)(ydots - 1));
595
596 /* bfxdel2 = (bfx3rd - bfxmin)/(ydots-1) */
597 sub_bf(bfxdel2, bfx3rd, bfxmin);
598 div_a_bf_int(bfxdel2, (U16)(ydots - 1));
599
600 /* bfydel2 = (bfy3rd - bfymin)/(xdots-1) */
601 sub_bf(bfydel2, bfy3rd, bfymin);
602 div_a_bf_int(bfydel2, (U16)(xdots - 1));
603
604 abs_bf(bfclosenuff,bfxdel);
605 if(cmp_bf(abs_bf(bftemp1,bfxdel2),bfclosenuff) > 0)
606 copy_bf(bfclosenuff,bftemp1);
607 if(cmp_bf(abs_bf(bftemp1,bfydel),abs_bf(bftemp2,bfydel2)) > 0)
608 {
609 if(cmp_bf(bftemp1,bfclosenuff) > 0)
610 copy_bf(bfclosenuff,bftemp1);
611 }
612 else if(cmp_bf(bftemp2,bfclosenuff) > 0)
613 copy_bf(bfclosenuff,bftemp2);
614 {
615 int t;
616 t = abs(periodicitycheck);
617 while(t--)
618 half_a_bf(bfclosenuff);
619 }
620
621 c_exp = (int)param[2];
622 switch (fractype)
623 {
624 case JULIAFP:
625 copy_bf(bfparm.x, bfparms[0]);
626 copy_bf(bfparm.y, bfparms[1]);
627 break;
628 case FPMANDELZPOWER:
629 init_big_pi();
630 if((double)c_exp == param[2] && (c_exp & 1)) /* odd exponents */
631 symmetry = XYAXIS_NOPARM;
632 if(param[3] != 0)
633 symmetry = NOSYM;
634 break;
635 case FPJULIAZPOWER:
636 init_big_pi();
637 copy_bf(bfparm.x, bfparms[0]);
638 copy_bf(bfparm.y, bfparms[1]);
639 if((c_exp & 1) || param[3] != 0.0 || (double)c_exp != param[2] )
640 symmetry = NOSYM;
641 break;
642 }
643
644 restore_stack(saved);
645 return (1);
646 }
647
648 #if (_MSC_VER >= 700)
649 #pragma code_seg ( ) /* back to normal segment */
650 #endif
651
652 int mandelbn_per_pixel()
653 {
654 /* parm.x = xxmin + col*delx + row*delx2 */
655 mult_bn_int(bnparm.x, bnxdel, (U16)col);
656 mult_bn_int(bntmp, bnxdel2, (U16)row);
657
658 add_a_bn(bnparm.x, bntmp);
659 add_a_bn(bnparm.x, bnxmin);
660
661 /* parm.y = yymax - row*dely - col*dely2; */
662 /* note: in next four lines, bnold is just used as a temporary variable */
663 mult_bn_int(bnold.x, bnydel, (U16)row);
664 mult_bn_int(bnold.y, bnydel2, (U16)col);
665 add_a_bn(bnold.x, bnold.y);
666 sub_bn(bnparm.y, bnymax, bnold.x);
667
668 copy_bn(bnold.x, bnparm.x);
669 copy_bn(bnold.y, bnparm.y);
670
671 if((inside == BOF60 || inside == BOF61) && !nobof)
672 {
673 /* kludge to match "Beauty of Fractals" picture since we start
674 Mandelbrot iteration with init rather than 0 */
675 floattobn(bnold.x,param[0]); /* initial pertubation of parameters set */
676 floattobn(bnold.y,param[1]);
677 coloriter = -1;
678 }
679 else
680 {
681 floattobn(bnnew.x,param[0]);
682 floattobn(bnnew.y,param[1]);
683 add_a_bn(bnold.x,bnnew.x);
684 add_a_bn(bnold.y,bnnew.y);
685 }
686
687 /* square has side effect - must copy first */
688 copy_bn(bnnew.x, bnold.x);
689 copy_bn(bnnew.y, bnold.y);
690
691 /* Square these to rlength bytes of precision */
692 square_bn(bntmpsqrx, bnnew.x);
693 square_bn(bntmpsqry, bnnew.y);
694
695 return (1); /* 1st iteration has been done */
696 }
697
698 int mandelbf_per_pixel()
699 {
700 /* parm.x = xxmin + col*delx + row*delx2 */
701 mult_bf_int(bfparm.x, bfxdel, (U16)col);
702 mult_bf_int(bftmp, bfxdel2, (U16)row);
703
704 add_a_bf(bfparm.x, bftmp);
705 add_a_bf(bfparm.x, bfxmin);
706
707 /* parm.y = yymax - row*dely - col*dely2; */
708 /* note: in next four lines, bfold is just used as a temporary variable */
709 mult_bf_int(bfold.x, bfydel, (U16)row);
710 mult_bf_int(bfold.y, bfydel2, (U16)col);
711 add_a_bf(bfold.x, bfold.y);
712 sub_bf(bfparm.y, bfymax, bfold.x);
713
714 copy_bf(bfold.x, bfparm.x);
715 copy_bf(bfold.y, bfparm.y);
716
717 if((inside == BOF60 || inside == BOF61) && !nobof)
718 {
719 /* kludge to match "Beauty of Fractals" picture since we start
720 Mandelbrot iteration with init rather than 0 */
721 floattobf(bfold.x,param[0]); /* initial pertubation of parameters set */
722 floattobf(bfold.y,param[1]);
723 coloriter = -1;
724 }
725 else
726 {
727 floattobf(bfnew.x,param[0]);
728 floattobf(bfnew.y,param[1]);
729 add_a_bf(bfold.x,bfnew.x);
730 add_a_bf(bfold.y,bfnew.y);
731 }
732
733 /* square has side effect - must copy first */
734 copy_bf(bfnew.x, bfold.x);
735 copy_bf(bfnew.y, bfold.y);
736
737 /* Square these to rbflength bytes of precision */
738 square_bf(bftmpsqrx, bfnew.x);
739 square_bf(bftmpsqry, bfnew.y);
740
741 return (1); /* 1st iteration has been done */
742 }
743
744 int
745 juliabn_per_pixel()
746 {
747 /* old.x = xxmin + col*delx + row*delx2 */
748 mult_bn_int(bnold.x, bnxdel, (U16)col);
749 mult_bn_int(bntmp, bnxdel2, (U16)row);
750
751 add_a_bn(bnold.x, bntmp);
752 add_a_bn(bnold.x, bnxmin);
753
754 /* old.y = yymax - row*dely - col*dely2; */
755 /* note: in next four lines, bnnew is just used as a temporary variable */
756 mult_bn_int(bnnew.x, bnydel, (U16)row);
757 mult_bn_int(bnnew.y, bnydel2, (U16)col);
758 add_a_bn(bnnew.x, bnnew.y);
759 sub_bn(bnold.y, bnymax, bnnew.x);
760
761 /* square has side effect - must copy first */
762 copy_bn(bnnew.x, bnold.x);
763 copy_bn(bnnew.y, bnold.y);
764
765 /* Square these to rlength bytes of precision */
766 square_bn(bntmpsqrx, bnnew.x);
767 square_bn(bntmpsqry, bnnew.y);
768
769 return (1); /* 1st iteration has been done */
770 }
771
772 int
773 juliabf_per_pixel()
774 {
775 /* old.x = xxmin + col*delx + row*delx2 */
776 mult_bf_int(bfold.x, bfxdel, (U16)col);
777 mult_bf_int(bftmp, bfxdel2, (U16)row);
778
779 add_a_bf(bfold.x, bftmp);
780 add_a_bf(bfold.x, bfxmin);
781
782 /* old.y = yymax - row*dely - col*dely2; */
783 /* note: in next four lines, bfnew is just used as a temporary variable */
784 mult_bf_int(bfnew.x, bfydel, (U16)row);
785 mult_bf_int(bfnew.y, bfydel2, (U16)col);
786 add_a_bf(bfnew.x, bfnew.y);
787 sub_bf(bfold.y, bfymax, bfnew.x);
788
789 /* square has side effect - must copy first */
790 copy_bf(bfnew.x, bfold.x);
791 copy_bf(bfnew.y, bfold.y);
792
793 /* Square these to rbflength bytes of precision */
794 square_bf(bftmpsqrx, bfnew.x);
795 square_bf(bftmpsqry, bfnew.y);
796
797 return (1); /* 1st iteration has been done */
798 }
799
800 int
801 JuliabnFractal()
802 {
803 /* Don't forget, with bn_t numbers, after multiplying or squaring */
804 /* you must shift over by shiftfactor to get the bn number. */
805
806 /* bntmpsqrx and bntmpsqry were previously squared before getting to */
807 /* this function, so they must be shifted. */
808
809 /* new.x = tmpsqrx - tmpsqry + parm.x; */
810 sub_a_bn(bntmpsqrx+shiftfactor, bntmpsqry+shiftfactor);
811 add_bn(bnnew.x, bntmpsqrx+shiftfactor, bnparm.x);
812
813 /* new.y = 2 * bnold.x * bnold.y + parm.y; */
814 mult_bn(bntmp, bnold.x, bnold.y); /* ok to use unsafe here */
815 double_a_bn(bntmp+shiftfactor);
816 add_bn(bnnew.y, bntmp+shiftfactor, bnparm.y);
817
818 return bignumbailout();
819 }
820
821 int
822 JuliabfFractal()
823 {
824 /* new.x = tmpsqrx - tmpsqry + parm.x; */
825 sub_a_bf(bftmpsqrx, bftmpsqry);
826 add_bf(bfnew.x, bftmpsqrx, bfparm.x);
827
828 /* new.y = 2 * bfold.x * bfold.y + parm.y; */
829 mult_bf(bftmp, bfold.x, bfold.y); /* ok to use unsafe here */
830 double_a_bf(bftmp);
831 add_bf(bfnew.y, bftmp, bfparm.y);
832 return bigfltbailout();
833 }
834
835 int
836 JuliaZpowerbnFractal()
837 {
838 _BNCMPLX parm2;
839 int saved; saved = save_stack();
840
841 parm2.x = alloc_stack(bnlength);
842 parm2.y = alloc_stack(bnlength);
843
844 floattobn(parm2.x,param[2]);
845 floattobn(parm2.y,param[3]);
846 ComplexPower_bn(&bnnew,&bnold,&parm2);
847 add_bn(bnnew.x,bnparm.x,bnnew.x+shiftfactor);
848 add_bn(bnnew.y,bnparm.y,bnnew.y+shiftfactor);
849 restore_stack(saved);
850 return bignumbailout();
851 }
852
853 int
854 JuliaZpowerbfFractal()
855 {
856 _BFCMPLX parm2;
857 int saved; saved = save_stack();
858
859 parm2.x = alloc_stack(bflength+2);
860 parm2.y = alloc_stack(bflength+2);
861
862 floattobf(parm2.x,param[2]);
863 floattobf(parm2.y,param[3]);
864 ComplexPower_bf(&bfnew,&bfold,&parm2);
865 add_bf(bfnew.x,bfparm.x,bfnew.x);
866 add_bf(bfnew.y,bfparm.y,bfnew.y);
867 restore_stack(saved);
868 return bigfltbailout();
869 }
870
871
872 #if 0
873 /*
874 the following is an example of how you can take advantage of the bn_t
875 format to squeeze a little more precision out of the calculations.
876 */
877 int
878 JuliabnFractal()
879 {
880 int oldbnlength;
881 bn_t mod;
882 /* using partial precision multiplications */
883
884 /* bnnew.x = bntmpsqrx - bntmpsqry + bnparm.x; */
885 /*
886 * Since tmpsqrx and tmpsqry where just calculated to rlength bytes of
887 * precision, we might as well keep that extra precision in this next
888 * subtraction. Therefore, use rlength as the length.
889 */
890
891 oldbnlength = bnlength;
892 bnlength = rlength; sub_a_bn(bntmpsqrx, bntmpsqry); bnlength = oldbnlength;
893
894 /*
895 * Now that bntmpsqry has been sutracted from bntmpsqrx, we need to treat
896 * tmpsqrx as a single width bignumber, so shift to bntmpsqrx+shiftfactor.
897 */
898 add_bn(bnnew.x, bntmpsqrx + shiftfactor, bnparm.x);
899
900 /* new.y = 2 * bnold.x * bnold.y + old.y; */
901 /* Multiply bnold.x*bnold.y to rlength precision. */
902 mult_bn(bntmp, bnold.x, bnold.y);
903
904 /*
905 * Double bnold.x*bnold.y by shifting bits, including one of those bits
906 * calculated in the previous mult_bn(). Therefore, use rlength.
907 */
908 bnlength = rlength; double_a_bn(bntmp); bnlength = oldbnlength;
909
910 /* Convert back to a single width bignumber and add bnparm.y */
911 add_bn(bnnew.y, bntmp + shiftfactor, bnparm.y);
912
913 copy_bn(bnold.x, bnnew.x);
914 copy_bn(bnold.y, bnnew.y);
915
916 /* Square these to rlength bytes of precision */
917 square_bn(bntmpsqrx, bnold.x);
918 square_bn(bntmpsqry, bnold.y);
919
920 /* And add the full rlength precision to get those extra bytes */
921 bnlength = rlength;
922 add_bn(bntmp, bntmpsqrx, bntmpsqry);
923 bnlength = oldbnlength;
924
925 mod = bntmp + (rlength) - (intlength << 1); /* where int part starts
926 * after mult */
927 /*
928 * equivalent to, but faster than, mod = bn_int(tmp+shiftfactor);
929 */
930
931 if ((magnitude = *mod) >= rqlim)
932 return (1);
933 return (0);
934 } 935 #endif
936
937 _CMPLX cmplxbntofloat(_BNCMPLX *s)
938 {
939 _CMPLX t;
940 t.x = (double)bntofloat(s->x);
941 t.y = (double)bntofloat(s->y);
942 return(t);
943 }
944
945 _CMPLX cmplxbftofloat(_BFCMPLX *s)
946 {
947 _CMPLX t;
948 t.x = (double)bftofloat(s->x);
949 t.y = (double)bftofloat(s->y);
950 return(t);
951 }
952
953 _BFCMPLX *cmplxlog_bf(_BFCMPLX *t, _BFCMPLX *s)
954 {
955 square_bf(t->x,s->x);
956 square_bf(t->y,s->y);
957 add_a_bf(t->x,t->y);
958 ln_bf(t->x,t->x);
959 half_a_bf(t->x);
960 atan2_bf(t->y,s->y,s->x);
961 return(t);
962 }
963
964 _BFCMPLX *cplxmul_bf( _BFCMPLX *t, _BFCMPLX *x, _BFCMPLX *y)
965 {
966 bf_t tmp1;
967 int saved; saved = save_stack();
968 tmp1 = alloc_stack(rbflength+2);
969 mult_bf(t->x, x->x, y->x);
970 mult_bf(t->y, x->y, y->y);
971 sub_bf(t->x,t->x,t->y);
972
973 mult_bf(tmp1, x->x, y->y);
974 mult_bf(t->y, x->y, y->x);
975 add_bf(t->y,tmp1,t->y);
976 restore_stack(saved);
977 return(t);
978 }
979
980 _BFCMPLX *ComplexPower_bf(_BFCMPLX *t, _BFCMPLX *xx, _BFCMPLX *yy)
981 {
982 _BFCMPLX tmp;
983 bf_t e2x, siny, cosy;
984 int saved; saved = save_stack();
985 e2x = alloc_stack(rbflength+2);
986 siny = alloc_stack(rbflength+2);
987 cosy = alloc_stack(rbflength+2);
988 tmp.x = alloc_stack(rbflength+2);
989 tmp.y = alloc_stack(rbflength+2);
990
991 /* 0 raised to anything is 0 */
992 if (is_bf_zero(xx->x) && is_bf_zero(xx->y))
993 {
994 clear_bf(t->x);
995 clear_bf(t->y);
996 return(t);
997 }
998
999 cmplxlog_bf(t, xx);
1000 cplxmul_bf(&tmp, t, yy);
1001 exp_bf(e2x,tmp.x);
1002 sincos_bf(siny,cosy,tmp.y);
1003 mult_bf(t->x, e2x, cosy);
1004 mult_bf(t->y, e2x, siny);
1005 restore_stack(saved);
1006 return(t);
1007 }
1008
1009 _BNCMPLX *cmplxlog_bn(_BNCMPLX *t, _BNCMPLX *s)
1010 {
1011 square_bn(t->x,s->x);
1012 square_bn(t->y,s->y);
1013 add_a_bn(t->x+shiftfactor,t->y+shiftfactor);
1014 ln_bn(t->x,t->x+shiftfactor);
1015 half_a_bn(t->x);
1016 atan2_bn(t->y,s->y,s->x);
1017 return(t);
1018 }
1019
1020 _BNCMPLX *cplxmul_bn( _BNCMPLX *t, _BNCMPLX *x, _BNCMPLX *y)
1021 {
1022 bn_t tmp1;
1023 int saved; saved = save_stack();
1024 tmp1 = alloc_stack(rlength);
1025 mult_bn(t->x, x->x, y->x);
1026 mult_bn(t->y, x->y, y->y);
1027 sub_bn(t->x,t->x+shiftfactor,t->y+shiftfactor);
1028
1029 mult_bn(tmp1, x->x, y->y);
1030 mult_bn(t->y, x->y, y->x);
1031 add_bn(t->y,tmp1+shiftfactor,t->y+shiftfactor);
1032 restore_stack(saved);
1033 return(t);
1034 }
1035
1036 /* note: ComplexPower_bn() returns need to be +shiftfactor'ed */
1037 _BNCMPLX *ComplexPower_bn(_BNCMPLX *t, _BNCMPLX *xx, _BNCMPLX *yy)
1038 {
1039 _BNCMPLX tmp;
1040 bn_t e2x, siny, cosy;
1041 int saved; saved = save_stack();
1042 e2x = alloc_stack(bnlength);
1043 siny = alloc_stack(bnlength);
1044 cosy = alloc_stack(bnlength);
1045 tmp.x = alloc_stack(rlength);
1046 tmp.y = alloc_stack(rlength);
1047
1048 /* 0 raised to anything is 0 */
1049 if (is_bn_zero(xx->x) && is_bn_zero(xx->y))
1050 {
1051 clear_bn(t->x);
1052 clear_bn(t->y);
1053 return(t);
1054 }
1055
1056 cmplxlog_bn(t, xx);
1057 cplxmul_bn(&tmp, t, yy);
1058 exp_bn(e2x,tmp.x);
1059 sincos_bn(siny,cosy,tmp.y);
1060 mult_bn(t->x, e2x, cosy);
1061 mult_bn(t->y, e2x, siny);
1062 restore_stack(saved);
1063 return(t);
1064 }
1065