Looping in the math evaluator?

SlavaCat118 · May 12, 2023, 3:27am

I’m looking to track how many times I must multiply x by y such that the product is a factor of x + y:

"n=1;(x*(y^n))%(x+y)==0"

I’ve looking into directly evaluating n, but it doesn’t seem too fun.
How should I record how many times it takes to increment n by 1 (starting from n=1) until the equation evaluates to true? I’d like it to be part of a single fill expression, but I’m open to more verbose solutions.

As a more minimal example:

"a=0;a==5"

How do I get the expression to return how many times it takes to increment a by 1 till a==5

Reptorian · May 12, 2023, 5:14am

I think you will end up in cases where the loop will never seemingly terminate. Like for x=5 and y=6. You’re going to need a mathematician to solve this. Are you sure you need to use ^? Edit: Tested and see no solution for this case though I think it might exist as a huge number.

It might help to know what this is about.

In the simple case, it seems to be:

n=0;
v=5;
repeat(inf,
 if(++n==v,break(););
);
count=n;

Note: Not tested.

Also, this might help - Modular Exponentiation (Power in Modular Arithmetic) - GeeksforGeeks

David_Tschumperle · May 12, 2023, 6:30am

Ahah, nice problem, I like it.

Some remarks:

if y==0 , your condition becomes 0%x == 0 which is true for all x (except undefined case when x==0).
if y==1, your condition becomes x%(x + 1)==0 which never holds except for x==0.

So, the cases y==0 and y==1 are not super interesting
For other y, I suggest this:

foo :
  300,300
  f "y<=1?0:( for (n = 1, (x*(y^n))%(x+y)!=0 && n<300, ++n); n)"

which gives:
res1

(not visually super interesting ).

Maybe with a log added ?

res2

a bit better, but still I wonder why you are looking for this kind of images?

snibgo · May 12, 2023, 7:34am

I’m looking to track how many times I must multiply x by y such that the product is a factor of x + y:

"n=1;(x*(y^n))%(x+y)==0"

My slightly-mathematician opinion: (EDIT: assuming x, y and n are integers) This is true if any only if x==0 or y==0. And then it is true for all values of n.

Reptorian · May 12, 2023, 5:06pm

I want to point out that this isn’t accurate because of precision and limited int of G’MIC. And I think there might be some more rules here and there.

I tested with Python which has “infinite int support”:

x=15
y=5
sum_xy=x+y

for n in range(1,1000+1):
    if ((x*(y**n))%sum_xy)==0:
        a=((x*(y**n)))
        print(((x*(y**n))))
        print(sum_xy)
        print(a//sum_xy)
        print(n)
        break
    
print("Finished")

Also, forgot to mention, regarding precision, see my rep_num2base_float_max code to see how one could get around the precision issue. This works around the limitation of int within G’MIC by splitting the number being represented as the maximum known number that can be represented in IEEE-754(?) as the base.

#@cli rep_num2base_float_max: number
#@cli : Create a image which represent a positive integer number in base 16777216. Used for emulating BIGINT process.
+rep_num2base_float_max:
user_number,base_float=$1,{1<<24}

('$1') fill[-1] inrange(i,_'0',_'9',1,1) if !im error inv_char_det fi rm.

1

for $user_number>9007199254740992 # The number seen here is equal to 1<<53

 eval "
  const _0=_'0';
  dividend='"$user_number"';
  res=0;
  repeat(size(dividend),p,
   res=(res*10+(dividend[p]-_0))%$base_float;
  );
  da_push(#-1,res);
  "

 1

 eval "
  number='"$user_number"';
  const divisor=$base_float;
  const n_size=size(number);
  const _0=_'0';
  index=0;
  dividend=number[index]-_0;

  while(dividend<divisor,
   dividend=dividend*10+(number[++index]-_0);
  );

  while(n_size>index,
   da_push(#-1,int(dividend/divisor)+_0);
   dividend=(dividend%divisor)*10+(number[++index]-_0);
  );

  da_freeze(#-1);
  "

 user_number={t}
 remove[-1]

done

eval "
 const base_float=$base_float;
 user_number=$user_number;

 do(
  da_push(#-1,user_number%base_float);
  user_number=int(user_number/base_float);
 ,user_number);

 resize(#-1,da_size(#-1),1,1,1,-1);"

David_Tschumperle · May 12, 2023, 5:18pm

You are right, G’MIC uses double-precision values, so when n becomes high, the modulo will incorrectly returns 0.
There must be something more fancy to do.

snibgo · May 12, 2023, 7:45pm

We don’t need to process large numbers.

Suppose we know the remainder after dividing A by C is R. In math-speak: A % C = R.

We want to know the remainder after dividing (A*y) by C. That remainder is (R*y)%C.

So to find each remainder, we don’t need to calculate (x*y^n). We need only to multiply the remainder by y, and the remainder is always in the range [0,x+y-1].

EDIT: A consequence of the above is that the remainder goes in cycles as we increase n. If one remainder is R1 and the next is R2, then R1 is always followed by R2.

The length of the cycle is in the range [1,x+y-1].

Suppose the remainder of A%C is zero. Then A is exactly divisible by C. The next remainder is (A*y)%C, but (A*y) must also be exactly divisible by C. So is zero remainder is always followed by another zero remainder. Hence it has a cycle length of one. Hence if we have a non-zero remainder for some x and y, we will never have a zero remainder, however many times we multiply by y.

Hence to the original question “how many times I must multiply x by y such that the product is a factor of x + y”, the answer is “Try it just once. If the remainder is zero or non-zero, it will remain so however many times we multiply by y.”

EDIT. No, this is wrong, somewhere. Sorry.

Reptorian · May 12, 2023, 9:19pm

Trying it once doesn’t work for in case of x=6, and y=10. It takes 3 cycles in this case.

Python foo example:

import math

x=6
y=10
sum_xy=x+y

for n in range(1,500+1):
    if ((x*(y**n))%sum_xy)==0:
        a=((x*(y**n)))
        print(f'x:{x}')
        print(f'y:{y}')
        print(f'xy_sum:{sum_xy}')
        print(f'((x*(y**n))):{((x*(y**n)))}')
        print(f'gcd:{math.gcd(x,y,sum_xy)}')   
        print(f'a/n:{a//sum_xy}')
        print(n)
        break
    
print("Finished")

Do you have a code to clarify on what you’re thinking of?

SlavaCat118 · May 12, 2023, 10:23pm

Thanks! I wasn’t aware of for loops in this context.

I was messing around, experimenting, and ran into the problem of needing to iterate a variable until a condition was true (I guess saying it like this makes the for loop obvious). The point isn’t really about the images, hence the name of the post and the absolutely minimal example, I was just trying to figure out this problem:

I included my application of it incase I was going about it the wrong way. I agree the images are rather boring, but hopefully with this knowledge I can make something slightly more interesting.

EDIT:
Something like this

300,300,1,1,"(x*y)==0?a=10:a=(x*y); (for(n=1, a%n==0 && n<100, ++n); n)"

Reptorian · May 12, 2023, 11:18pm

If you’re looking into experimenting with loops, you can look into knowing 4 different loops.

You learned for loop.

Here are other loops within math parser:

do(,while_condition);
repeat(n,index_name,); ( can also be repeat(n,); )
while(condition,);

repeat(n…); is the most common of them all because most of the time, n can be known.

Outside the math parser, for is more limited as in it’s just a condition. There is do … while …, and repeat … done.

snibgo · May 13, 2023, 10:23am

Sorry for the thread drift. Here is a python script for the original problem. If I have the maths and python correct, it will always either find a remainder that has occurred before (and therefore we are in a cycle), or a remainder of zero (and hence the next remainder would be zero). It should never say “Not finished”.

The script finds the modulus of remainders, so doesn’t need large numbers. It uses a python set to keep track of what remainders have been found.

x=5
y=6

rems=set(())

sum_xy=x+y

remainder=(x*y)%sum_xy
print("remainder =", remainder)
rems.add(remainder)

finished=False

for n in range(1,sum_xy+1):
    remainder=(remainder*y)%sum_xy
    print("remainder =", remainder)
    if remainder==0:
        print("Found remainder zero at n=", n)
        finished=True
        break
    if remainder in rems:
        print("Already found remainder ", remainder, "at n=", n)
        finished=True
        break
    rems.add(remainder)
    
print("Done. rems = ", rems)

if not finished:
    print("Not finished")

David_Tschumperle · May 13, 2023, 1:57pm

If I understand @snibgo correctly, this could be done like this in G’MIC :

foo :
  300,300
  f "const itermax = 100;
     S = x + y;
     P = x;
     S?repeat (itermax,n,
       P = (P*y)%S;
       !P?break();
     ):0;
     n>=itermax?0:(n+1)"

which returns this image (here, normalized in [0,255], but actual range is [0,8]):

Is this correct?

David_Tschumperle · May 13, 2023, 2:06pm

Now the fun starts:
If I compute this for a 4000x4000 image, and just output each pixel to 0 or 1 (1 meaning : there exists an integer n s.t. (x*y^n)%(x+y)==0, then the image looks interesting !
(I’ve reduced it to 1000x1000 afterwards).

Reptorian · May 13, 2023, 2:38pm

I understood as more like this:

$ 300,300 rep_slavacat_gen_art
rep_slavacat_gen_art:
if !$! error "n_imgs>0==F" fi

1x$_cpus

fill[-{$_cpus+1}] :"
 img_pos=k+t+1;
 sum_xy=x+y;
 remainder=(x*y)%sum_xy;
 da_push(#img_pos,remainder);
 finished=0;
 repeat(sum_xy,n,
  remainder=(remainder*y)%sum_xy;
  if(!remainder,
   da_remove(#img_pos,0,da_size(#img_pos)-1);
   finished=1;
   break();
  );
  repeat(da_size(#img_pos),p,
   if(remainder==i[#img_pos,p],
    finished=1;
    da_remove(#img_pos,0,da_size(#img_pos)-1);
    break();
   );
  );
  da_push(#img_pos,remainder);
  if(finished,break(););
 );
 if(!finished,
  da_remove(#img_pos,0,da_size(#img_pos)-1);
 );
 n;
 img_pos;
 "

It seem to be a more faithful code given the similarities.

It does not work unfortunately. Hence my last line which is debugging.

EDIT:

Considering multi-threaded processing using multiple images as dynamic array is bugged, I resorted to single thread. Here’s the result:

The code:

rep_slavacat_gen_art:
if !$! error "n_imgs>0==F" fi

1

fill[-2] >"
 sum_xy=x+y;
 remainder=(x*y)%max(1,sum_xy);
 da_push(#-1,remainder);
 finished=0;
 repeat(sum_xy,n,
  remainder=(remainder*y)%sum_xy;
  if(!remainder,
   da_remove(#-1,0,da_size(#-1)-1);
   finished=1;
   break();
  );
  repeat(da_size(#-1),p,
   if(remainder==i[#-1,p],
    finished=1;
    da_remove(#-1,0,da_size(#-1)-1);
    break();
   );
  );
  da_push(#-1,remainder);
  if(finished,
   da_remove(#-1,0,da_size(#-1)-1);
   break();
  );
 );
 if(!finished,
  da_remove(#-1,0,da_size(#-1)-1);
 );
 n;
 "

Comparing David’s result with my own Python code, it is determined David’s code works.

snibgo · May 13, 2023, 8:10pm

Yes, that’s it. Doing it in ImageMagick, as a bash command:

magick \
  -size 300x300 xc:Black \
  -fx " itermax=100;
        sumxy = i + j;
        remainder = i;
        for ( nn = 1,
              nn < itermax && remainder != 0,
              nn++; remainder = (remainder * j) % sumxy
            );
        nn >= itermax ? 0 : nn+1
      " \
  -auto-level \
  x.png

This is black or non-black where David’s is the same.

Reptorian · May 13, 2023, 8:51pm

One last thing, I found out that you don’t need to iterate 100 times.

The maximum number of iteration you need is equal to max(2,int(log2(max(x,y)))) or even add that with 1 for good measure as I have no proof that this is true, but it seems to be the case. At the moment, it’s only a conjecture. You can even confirm it by looking at iM in G’MIC and see your image size.

snibgo · May 13, 2023, 9:59pm

That would be useful, if it is true. It doesn’t seem to be true (or I misunderstand you).

For the case of x=5, y=6:

remainder = 8
remainder = 4
remainder = 2
remainder = 1
remainder = 6
remainder = 3
remainder = 7
remainder = 9
remainder = 10
remainder = 5
remainder = 8
Already found remainder  8 at n= 10
Done. rems =  {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}

The cycle length is 10. But max(2,int(log2(max(x,y)))) = 2.

There is an upper limit of (x+y). That number of iterations will always be enough to either find a zero or a repeated remainder. For example, with x=5 and y=6, 11 iterations are certain to be enough. It turns out that we need 10 (because zero doesn’t occur).

Reptorian · May 13, 2023, 11:35pm

You misunderstand me, I was using David Tschumperle’s code, and this is the test for the conjecture:

rep_slavacat_gen_art:
if !$! error "n_imgs>0==F" fi

+f[-1] "const itermax = 100;
 S = x + y;
 P = x;
 S?repeat (itermax,n,
  P = (P*y)%S;
  !P?break();
 ):0;
 n>=itermax?0:(n+1)"

+f[-2] "itermax = (r=max(x,y);r?max(2,int(log2(r))));
 S = x + y;
 P = x;
 S?repeat(itermax,n,
  P = (P*y)%S;
  !P?break();
 ):0;
 n>=itermax?0:(n+1)"
 
+blend[-2,-1] difference gt. 0

If you execute this cli code:

$ size=4098 $size,$size rep_slavacat_gen_art echo {int(log2($size))} if !iM echo Pass! fi

Then, you’ll notice it keeps passing no matter how much you change size. So, I can only conclude my conjecture is correct in case of your ImageMagick code as a bash command or David’s command.

There is also a massive speed of difference compared with const itermax=100;

C:\Windows\System32>gmic size=8192 $size,$size tic rep_slavacat_gen_art toc
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable 'size=8192'.
[gmic]-0./ Input black image at position 0 (1 image 8192x8192x1x1).
[gmic]-1./ Initialize timer.
[gmic]-1./ Elapsed time: 3.266 s.
[gmic]-1./ Display image [0] = '[unnamed]'.
[0] = '[unnamed]':
  size = (8192,8192,1,1) [256 Mio of float32].
  data = (0,1,1,1,1,1,1,1,1,1,1,1,(...),0,0,0,0,0,0,0,0,0,0,0,0).
  min = 0, max = 12, mean = 0.00572254, std = 0.171757, coords_min = (0,0,0,0), coords_max = (8190,2,0,0).
[gmic]-1./ End G'MIC interpreter.

C:\Windows\System32>gmic size=8192 $size,$size tic rep_slavacat_gen_art toc
[gmic]-0./ Start G'MIC interpreter.
[gmic]-0./ Set local variable 'size=8192'.
[gmic]-0./ Input black image at position 0 (1 image 8192x8192x1x1).
[gmic]-1./ Initialize timer.
[gmic]-1./ Elapsed time: 22.926 s.
[gmic]-1./ Display image [0] = '[unnamed]'.
[0] = '[unnamed]':
  size = (8192,8192,1,1) [256 Mio of float32].
  data = (0,1,1,1,1,1,1,1,1,1,1,1,(...),0,0,0,0,0,0,0,0,0,0,0,0).
  min = 0, max = 12, mean = 0.00572254, std = 0.171757, coords_min = (0,0,0,0), coords_max = (8190,2,0,0).
[gmic]-1./ End G'MIC interpreter.