 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
DanRRight
Joined: 10 Mar 2008
Posts: 2863
Location: South Pole, Antarctica
|
 Posted: Sun Sep 15, 2024 4:11 am Post subject: Any MATLAB Pros are here? |
 |
|
Can you help to convert this Fortran code to Matlab?
| Code: | integer,parameter::k=selected_int_kind(18)
integer(k) v
integer U,L
111 READ(*,*) v
U = ishft(v, -24)
L = iand(v, Z'000000ffffff')
print*,L,U
! getting v back
v = int(L,k) + ishft(int(U,k),24)
print*,v
goto 111
end |
We discussed this code a year ago, see https://forums.silverfrost.com/viewtopic.php?t=4840.
This code reads from screen 64bit number v you put from the keyboard and splits it into two 32bit single precision numbers U and L. Then it combines two numbers into v back to check if all went OK.
This code takes only first 24bits of the v number to restrict U and L to values below approximately 20 million, because if you convert more and then assign these U and L integer numbers to real single precision 32bit numbers (and we have to do that), all further bits will be lost (mantissa of real number is just 24 bits or so). This code with this tricky way can split v numbers up to 10^15 instead of just 20 million or so |
|
| Back to top |
|
 |
DanRRight
Joined: 10 Mar 2008
Posts: 2863
Location: South Pole, Antarctica
|
| Posted: Mon Sep 16, 2024 3:45 pm Post subject: |
|
|
| Seems there was not a Matlab but the FTN95 problem. We wrote Matlab code but it produced different result for U. I just tried Gfortran and also one additional Fortran compiler and both also "agree to disagree" with FTN95. If we try any v number less than 16 million FTN95 gives U=1 but others having U=0. All generate the same correct value L. Fun is that FTN95 somehow still recovers v value back correctly, which means it has double problem ? Or it is all other compilers which have a problem? |
|
| Back to top |
|
|
PaulLaidler
Site Admin
Joined: 21 Feb 2005
Posts: 8011
Location: Salford, UK
|
| Posted: Mon Sep 16, 2024 4:16 pm Post subject: |
|
|
I don't have time to look at this now but I seem to recall that the last time we looked at this was http://forums.silverfrost.com/viewtopic.php?t=4840&highlight=ishft.
The Z constants required a kind specifier. There may be differences in the required Z form between one compiler and another. |
|
| Back to top |
|
|
DanRRight
Joined: 10 Mar 2008
Posts: 2863
Location: South Pole, Antarctica
|
| Posted: Tue Sep 17, 2024 6:22 am Post subject: |
|
|
That was good advice, thanks.
By the way Matlab simply refused to compile in this case giving an error. And Gfortran did not complain but gave correct result without me specifying 64bit Z number not as simply Znumber as Znumber_k.
May be it is still worth for FTN95 to make in this specific place an error because this trick is so obscure that only 1 per million people will find the reason of discrepancy. Problem also was that combining back the splitted numbers to check if splitting ended fine, was giving correct value - double error with correct income. Pure devilry |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2580
Location: Sydney
|
| Posted: Wed Sep 18, 2024 6:05 am Post subject: |
|
|
Dan,
I think the only problem is with iand, so I included a corrected I*8 mask.
Let me know if I have misunderstood.
I modified your code, to enforce a 64-bit integer mask, "mask24".
I tested for a number of values "v", but all worked with your 32-bit mask.
You didn't provide values for v. Do you have any values where v fails ?
Actually, any values larger than 2^55 will fail (U is -ve), but any others ?
Why use 24 bit split ?
This is my expanded code to identify any problems is:
| Code: | ! This works for all values smaller that 2^(24+31),
! ie no overflow on I84 U
integer, parameter :: i8=selected_int_kind(18)
integer(i8), parameter :: two = 2
integer(i8), parameter :: mask24 = two**24-1
integer(i8) :: vlist(6) = [ 16777216, 16777215, 134217728, 134217727, 1234567890, -1 ]
integer(i8) :: v
integer :: U,L, k
write (*,11) mask24, mask24, mask24
11 format ( 'Mask ',z20,i20,b33)
! vlist(5) = two ** 34 - 3 ! larger value
vlist(5) = two ** 58 - 3 ! larger value that will fail as longer than 55 bits
!111 READ(*,*) v
do k = 1,size(vlist)
! Old approach
v = vlist(k)
U = ishft ( v, -24 )
L = iand ( v, Z'000000ffffff' )
! print*,L,U
! getting v back
v = int(L,i8) + ishft(int(U,i8),24)
! print*,v
print*,'old',U,L,v,vlist(k)
! New approach with I*8 : mask24
v = vlist(k)
U = ishft ( v, -24 ) ! problem if U -ve
L = iand ( v, mask24 )
! print*,L,U
! getting v back
v = int(L,i8) + ishft(int(U,i8),24)
! print*,v
if ( v == vlist(k) ) then
print*,'new',U,L,v,vlist(k)
else
print*,'new',U,L,v,vlist(k),' ERROR'
end if
end do
end |
|
|
| Back to top |
|
|
DanRRight
Joined: 10 Mar 2008
Posts: 2863
Location: South Pole, Antarctica
|
| Posted: Sun Sep 29, 2024 7:37 pm Post subject: |
|
|
| 24 bit split used because after V splitted to two integer numbers U and L they will be assigned to real*4 numbers RU and RL and saved to media (that was done for speed considerations, RU and RL will be a part of some huge multi-GB or TB size arrays) . When file will be loaded again these RU and RL will be converted to V preserving all its first digits up to something like 10^15 value. This is not full Integer*8 but I hope that the 640...sorry I mean 10^15 64bit numbers will fit to everyone. How many Petabytes is it? |
|
| Back to top |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
