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